C ------------- BELOW IS DSYEVX -------------------- CAT > DSYEVX.F <<'CUT HERE............' SUBROUTINE DSYEVX( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU, $ ABSTOL, M, W, Z, LDZ, WORK, LWORK, IWORK, $ IFAIL, INFO ) * * -- LAPACK DRIVER ROUTINE (VERSION 1.1) -- * UNIV. OF TENNESSEE, UNIV. OF CALIFORNIA BERKELEY, NAG LTD., * COURANT INSTITUTE, ARGONNE NATIONAL LAB, AND RICE UNIVERSITY * MARCH 31, 1993 * * .. SCALAR ARGUMENTS .. CHARACTER JOBZ, RANGE, UPLO INTEGER IL, INFO, IU, LDA, LDZ, LWORK, M, N DOUBLE PRECISION ABSTOL, VL, VU * .. * .. ARRAY ARGUMENTS .. INTEGER IFAIL( * ), IWORK( * ) DOUBLE PRECISION A( LDA, * ), W( * ), WORK( * ), Z( LDZ, * ) * .. * * PURPOSE * ======= * * DSYEVX COMPUTES SELECTED EIGENVALUES AND, OPTIONALLY, EIGENVECTORS * OF A REAL SYMMETRIC MATRIX A. EIGENVALUES AND EIGENVECTORS CAN BE * SELECTED BY SPECIFYING EITHER A RANGE OF VALUES OR A RANGE OF INDICES * FOR THE DESIRED EIGENVALUES. * * ARGUMENTS * ========= * * JOBZ (INPUT) CHARACTER*1 * = 'N': COMPUTE EIGENVALUES ONLY; * = 'V': COMPUTE EIGENVALUES AND EIGENVECTORS. * * RANGE (INPUT) CHARACTER*1 * = 'A': ALL EIGENVALUES WILL BE FOUND. * = 'V': ALL EIGENVALUES IN THE HALF-OPEN INTERVAL (VL,VU] * WILL BE FOUND. * = 'I': THE IL-TH THROUGH IU-TH EIGENVALUES WILL BE FOUND. * * UPLO (INPUT) CHARACTER*1 * = 'U': UPPER TRIANGLE OF A IS STORED; * = 'L': LOWER TRIANGLE OF A IS STORED. * * N (INPUT) INTEGER * THE ORDER OF THE MATRIX A. N >= 0. * * A (INPUT/WORKSPACE) DOUBLE PRECISION ARRAY, DIMENSION (LDA, N) * ON ENTRY, THE SYMMETRIC MATRIX A. IF UPLO = 'U', THE * LEADING N-BY-N UPPER TRIANGULAR PART OF A CONTAINS THE * UPPER TRIANGULAR PART OF THE MATRIX A. IF UPLO = 'L', * THE LEADING N-BY-N LOWER TRIANGULAR PART OF A CONTAINS * THE LOWER TRIANGULAR PART OF THE MATRIX A. * ON EXIT, THE LOWER TRIANGLE (IF UPLO='L') OR THE UPPER * TRIANGLE (IF UPLO='U') OF A, INCLUDING THE DIAGONAL, IS * DESTROYED. * * LDA (INPUT) INTEGER * THE LEADING DIMENSION OF THE ARRAY A. LDA >= MAX(1,N). * * VL (INPUT) DOUBLE PRECISION * IF RANGE='V', THE LOWER BOUND OF THE INTERVAL TO BE SEARCHED * FOR EIGENVALUES. NOT REFERENCED IF RANGE = 'A' OR 'I'. * * VU (INPUT) DOUBLE PRECISION * IF RANGE='V', THE UPPER BOUND OF THE INTERVAL TO BE SEARCHED * FOR EIGENVALUES. NOT REFERENCED IF RANGE = 'A' OR 'I'. * * IL (INPUT) INTEGER * IF RANGE='I', THE INDEX (FROM SMALLEST TO LARGEST) OF THE * SMALLEST EIGENVALUE TO BE RETURNED. IL >= 1. * NOT REFERENCED IF RANGE = 'A' OR 'V'. * * IU (INPUT) INTEGER * IF RANGE='I', THE INDEX (FROM SMALLEST TO LARGEST) OF THE * LARGEST EIGENVALUE TO BE RETURNED. MIN(IL,N) <= IU <= N. * NOT REFERENCED IF RANGE = 'A' OR 'V'. * * ABSTOL (INPUT) DOUBLE PRECISION * THE ABSOLUTE ERROR TOLERANCE FOR THE EIGENVALUES. * AN APPROXIMATE EIGENVALUE IS ACCEPTED AS CONVERGED * WHEN IT IS DETERMINED TO LIE IN AN INTERVAL [A,B] * OF WIDTH LESS THAN OR EQUAL TO * * ABSTOL + EPS * MAX( |A|,|B| ) , * * WHERE EPS IS THE MACHINE PRECISION. IF ABSTOL IS LESS THAN * OR EQUAL TO ZERO, THEN EPS*|T| WILL BE USED IN ITS PLACE, * WHERE |T| IS THE 1-NORM OF THE TRIDIAGONAL MATRIX OBTAINED * BY REDUCING A TO TRIDIAGONAL FORM. * * SEE "COMPUTING SMALL SINGULAR VALUES OF BIDIAGONAL MATRICES * WITH GUARANTEED HIGH RELATIVE ACCURACY," BY DEMMEL AND * KAHAN, LAPACK WORKING NOTE #3. * * M (OUTPUT) INTEGER * THE TOTAL NUMBER OF EIGENVALUES FOUND. 0 <= M <= N. * IF RANGE = 'A', M = N, AND IF RANGE = 'I', M = IU-IL+1. * * W (OUTPUT) DOUBLE PRECISION ARRAY, DIMENSION (N) * ON NORMAL EXIT, THE FIRST M ENTRIES CONTAIN THE SELECTED * EIGENVALUES IN ASCENDING ORDER. * * Z (OUTPUT) DOUBLE PRECISION ARRAY, DIMENSION (LDZ, MAX(1,M)) * IF JOBZ = 'V', THEN IF INFO = 0, THE FIRST M COLUMNS OF Z * CONTAIN THE ORTHONORMAL EIGENVECTORS OF THE MATRIX * CORRESPONDING TO THE SELECTED EIGENVALUES. IF AN EIGENVECTOR * FAILS TO CONVERGE, THEN THAT COLUMN OF Z CONTAINS THE LATEST * APPROXIMATION TO THE EIGENVECTOR, AND THE INDEX OF THE * EIGENVECTOR IS RETURNED IN IFAIL. * IF JOBZ = 'N', THEN Z IS NOT REFERENCED. * NOTE: THE USER MUST ENSURE THAT AT LEAST MAX(1,M) COLUMNS ARE * SUPPLIED IN THE ARRAY Z; IF RANGE = 'V', THE EXACT VALUE OF M * IS NOT KNOWN IN ADVANCE AND AN UPPER BOUND MUST BE USED. * * LDZ (INPUT) INTEGER * THE LEADING DIMENSION OF THE ARRAY Z. LDZ >= 1, AND IF * JOBZ = 'V', LDZ >= MAX(1,N). * * WORK (WORKSPACE) DOUBLE PRECISION ARRAY, DIMENSION (LWORK) * ON EXIT, IF INFO = 0, WORK(1) RETURNS THE OPTIMAL LWORK. * * LWORK (INPUT) INTEGER * THE LENGTH OF THE ARRAY WORK. LWORK >= MAX(1,8*N). * FOR OPTIMAL EFFICIENCY, LWORK >= (NB+3)*N, * WHERE NB IS THE BLOCKSIZE FOR DSYTRD RETURNED BY ILAENV. * * IWORK (WORKSPACE) INTEGER ARRAY, DIMENSION (5*N) * * IFAIL (OUTPUT) INTEGER ARRAY, DIMENSION (N) * IF JOBZ = 'V', THEN IF INFO = 0, THE FIRST M ELEMENTS OF * IFAIL ARE ZERO. IF INFO > 0, THEN IFAIL CONTAINS THE * INDICES OF THE EIGENVECTORS THAT FAILED TO CONVERGE. * IF JOBZ = 'N', THEN IFAIL IS NOT REFERENCED. * * INFO (OUTPUT) INTEGER * = 0: SUCCESSFUL EXIT * < 0: IF INFO = -I, THE I-TH ARGUMENT HAD AN ILLEGAL VALUE * > 0: IF INFO = I, THEN I EIGENVECTORS FAILED TO CONVERGE. * THEIR INDICES ARE STORED IN ARRAY IFAIL. * * ===================================================================== * * .. PARAMETERS .. DOUBLE PRECISION ZERO, ONE PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) * .. * .. LOCAL SCALARS .. LOGICAL ALLEIG, INDEIG, LOWER, VALEIG, WANTZ CHARACTER ORDER INTEGER I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL, $ INDISP, INDIWO, INDTAU, INDWKN, INDWRK, ISCALE, $ ITMP1, J, JJ, LLWORK, LLWRKN, LOPT, NSPLIT DOUBLE PRECISION ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, $ SIGMA, SMLNUM, TMP1, VLL, VUU * .. * .. EXTERNAL FUNCTIONS .. LOGICAL LSAME DOUBLE PRECISION DLAMCH, DLANSY EXTERNAL LSAME, DLAMCH, DLANSY * .. * .. EXTERNAL SUBROUTINES .. EXTERNAL DCOPY, DLACPY, DORGTR, DORMTR, DSCAL, DSTEBZ, $ DSTEIN, DSTEQR, DSTERF, DSWAP, DSYTRD, XERBLA * .. * .. INTRINSIC FUNCTIONS .. INTRINSIC MAX, MIN, SQRT * .. * .. EXECUTABLE STATEMENTS .. * * TEST THE INPUT PARAMETERS. * LOWER = LSAME( UPLO, 'L' ) WANTZ = LSAME( JOBZ, 'V' ) ALLEIG = LSAME( RANGE, 'A' ) VALEIG = LSAME( RANGE, 'V' ) INDEIG = LSAME( RANGE, 'I' ) * INFO = 0 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN INFO = -1 ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN INFO = -2 ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN INFO = -3 ELSE IF( N.LT.0 ) THEN INFO = -4 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN INFO = -6 ELSE IF( VALEIG .AND. N.GT.0 .AND. VU.LE.VL ) THEN INFO = -8 ELSE IF( INDEIG .AND. IL.LT.1 ) THEN INFO = -9 ELSE IF( INDEIG .AND. ( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) ) THEN INFO = -10 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN INFO = -15 ELSE IF( LWORK.LT.MAX( 1, 8*N ) ) THEN INFO = -17 END IF * IF( INFO.NE.0 ) THEN CALL XERBLA( 'DSYEVX', -INFO ) RETURN END IF * * QUICK RETURN IF POSSIBLE * M = 0 IF( N.EQ.0 ) THEN WORK( 1 ) = 1 RETURN END IF * IF( N.EQ.1 ) THEN WORK( 1 ) = 7 IF( ALLEIG .OR. INDEIG ) THEN M = 1 W( 1 ) = A( 1, 1 ) ELSE IF( VL.LT.A( 1, 1 ) .AND. VU.GE.A( 1, 1 ) ) THEN M = 1 W( 1 ) = A( 1, 1 ) END IF END IF IF( WANTZ ) $ Z( 1, 1 ) = ONE RETURN END IF * * GET MACHINE CONSTANTS. * SAFMIN = DLAMCH( 'SAFE MINIMUM' ) EPS = DLAMCH( 'PRECISION' ) SMLNUM = SAFMIN / EPS BIGNUM = ONE / SMLNUM RMIN = SQRT( SMLNUM ) RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) ) * * SCALE MATRIX TO ALLOWABLE RANGE, IF NECESSARY. * ISCALE = 0 ABSTLL = ABSTOL IF( VALEIG ) THEN VLL = VL VUU = VU END IF ANRM = DLANSY( 'M', UPLO, N, A, LDA, WORK ) IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN ISCALE = 1 SIGMA = RMIN / ANRM ELSE IF( ANRM.GT.RMAX ) THEN ISCALE = 1 SIGMA = RMAX / ANRM END IF IF( ISCALE.EQ.1 ) THEN IF( LOWER ) THEN DO 10 J = 1, N CALL DSCAL( N-J+1, SIGMA, A( J, J ), 1 ) 10 CONTINUE ELSE DO 20 J = 1, N CALL DSCAL( J, SIGMA, A( 1, J ), 1 ) 20 CONTINUE END IF IF( ABSTOL.GT.0 ) $ ABSTLL = ABSTOL*SIGMA IF( VALEIG ) THEN VLL = VL*SIGMA VUU = VU*SIGMA END IF END IF * * CALL DSYTRD TO REDUCE SYMMETRIC MATRIX TO TRIDIAGONAL FORM. * INDTAU = 1 INDE = INDTAU + N INDD = INDE + N INDWRK = INDD + N LLWORK = LWORK - INDWRK + 1 CALL DSYTRD( UPLO, N, A, LDA, WORK( INDD ), WORK( INDE ), $ WORK( INDTAU ), WORK( INDWRK ), LLWORK, IINFO ) LOPT = 3*N + WORK( INDWRK ) * * IF ALL EIGENVALUES ARE DESIRED AND ABSTOL IS LESS THAN OR EQUAL TO * ZERO, THEN CALL DSTERF OR DORGTR AND SSTEQR. IF THIS FAILS FOR * SOME EIGENVALUE, THEN TRY DSTEBZ. * IF( ( ALLEIG .OR. ( INDEIG .AND. IL.EQ.1 .AND. IU.EQ.N ) ) .AND. $ ( ABSTOL.LE.ZERO ) ) THEN CALL DCOPY( N, WORK( INDD ), 1, W, 1 ) INDEE = INDWRK + 2*N IF( .NOT.WANTZ ) THEN CALL DCOPY( N-1, WORK( INDE ), 1, WORK( INDEE ), 1 ) CALL DSTERF( N, W, WORK( INDEE ), INFO ) ELSE CALL DLACPY( 'A', N, N, A, LDA, Z, LDZ ) CALL DORGTR( UPLO, N, Z, LDZ, WORK( INDTAU ), $ WORK( INDWRK ), LLWORK, IINFO ) CALL DCOPY( N-1, WORK( INDE ), 1, WORK( INDEE ), 1 ) CALL DSTEQR( JOBZ, N, W, WORK( INDEE ), Z, LDZ, $ WORK( INDWRK ), INFO ) IF( INFO.EQ.0 ) THEN DO 30 I = 1, N IFAIL( I ) = 0 30 CONTINUE END IF END IF IF( INFO.EQ.0 ) THEN M = N GO TO 40 END IF INFO = 0 END IF * * OTHERWISE, CALL DSTEBZ AND, IF EIGENVECTORS ARE DESIRED, SSTEIN. * IF( WANTZ ) THEN ORDER = 'B' ELSE ORDER = 'E' END IF INDIBL = 1 INDISP = INDIBL + N INDIWO = INDISP + N CALL DSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTLL, $ WORK( INDD ), WORK( INDE ), M, NSPLIT, W, $ IWORK( INDIBL ), IWORK( INDISP ), WORK( INDWRK ), $ IWORK( INDIWO ), INFO ) * IF( WANTZ ) THEN CALL DSTEIN( N, WORK( INDD ), WORK( INDE ), M, W, $ IWORK( INDIBL ), IWORK( INDISP ), Z, LDZ, $ WORK( INDWRK ), IWORK( INDIWO ), IFAIL, INFO ) * * APPLY ORTHOGONAL MATRIX USED IN REDUCTION TO TRIDIAGONAL * FORM TO EIGENVECTORS RETURNED BY DSTEIN. * INDWKN = INDE LLWRKN = LWORK - INDWKN + 1 CALL DORMTR( 'L', UPLO, 'N', N, M, A, LDA, WORK( INDTAU ), Z, $ LDZ, WORK( INDWKN ), LLWRKN, IINFO ) END IF * * IF MATRIX WAS SCALED, THEN RESCALE EIGENVALUES APPROPRIATELY. * 40 CONTINUE IF( ISCALE.EQ.1 ) THEN IF( INFO.EQ.0 ) THEN IMAX = M ELSE IMAX = INFO - 1 END IF CALL DSCAL( IMAX, ONE / SIGMA, W, 1 ) END IF * * IF EIGENVALUES ARE NOT IN ORDER, THEN SORT THEM, ALONG WITH * EIGENVECTORS. * IF( WANTZ ) THEN DO 60 J = 1, M - 1 I = 0 TMP1 = W( J ) DO 50 JJ = J + 1, M IF( W( JJ ).LT.TMP1 ) THEN I = JJ TMP1 = W( JJ ) END IF 50 CONTINUE * IF( I.NE.0 ) THEN ITMP1 = IWORK( INDIBL+I-1 ) W( I ) = W( J ) IWORK( INDIBL+I-1 ) = IWORK( INDIBL+J-1 ) W( J ) = TMP1 IWORK( INDIBL+J-1 ) = ITMP1 CALL DSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 ) IF( INFO.NE.0 ) THEN ITMP1 = IFAIL( I ) IFAIL( I ) = IFAIL( J ) IFAIL( J ) = ITMP1 END IF END IF 60 CONTINUE END IF * * SET WORK(1) TO OPTIMAL WORKSPACE SIZE. * WORK( 1 ) = MAX( 7*N, LOPT ) * RETURN * * END OF DSYEVX * END CUT HERE............ CAT > DORMTR.F <<'CUT HERE............' SUBROUTINE DORMTR( SIDE, UPLO, TRANS, M, N, A, LDA, TAU, C, LDC, $ WORK, LWORK, INFO ) * * -- LAPACK ROUTINE (VERSION 1.1) -- * UNIV. OF TENNESSEE, UNIV. OF CALIFORNIA BERKELEY, NAG LTD., * COURANT INSTITUTE, ARGONNE NATIONAL LAB, AND RICE UNIVERSITY * MARCH 31, 1993 * * .. SCALAR ARGUMENTS .. CHARACTER SIDE, TRANS, UPLO INTEGER INFO, LDA, LDC, LWORK, M, N * .. * .. ARRAY ARGUMENTS .. DOUBLE PRECISION A( LDA, * ), C( LDC, * ), TAU( * ), $ WORK( LWORK ) * .. * * PURPOSE * ======= * * DORMTR OVERWRITES THE GENERAL REAL M-BY-N MATRIX C WITH * * SIDE = 'L' SIDE = 'R' * TRANS = 'N': Q * C C * Q * TRANS = 'T': Q**T * C C * Q**T * * WHERE Q IS A REAL ORTHOGONAL MATRIX OF ORDER NQ, WITH NQ = M IF * SIDE = 'L' AND NQ = N IF SIDE = 'R'. Q IS DEFINED AS THE PRODUCT OF * NQ-1 ELEMENTARY REFLECTORS, AS RETURNED BY DSYTRD: * * IF UPLO = 'U', Q = H(NQ-1) . . . H(2) H(1); * * IF UPLO = 'L', Q = H(1) H(2) . . . H(NQ-1). * * ARGUMENTS * ========= * * SIDE (INPUT) CHARACTER*1 * = 'L': APPLY Q OR Q**T FROM THE LEFT; * = 'R': APPLY Q OR Q**T FROM THE RIGHT. * * UPLO (INPUT) CHARACTER*1 * = 'U': UPPER TRIANGLE OF A CONTAINS ELEMENTARY REFLECTORS * FROM DSYTRD; * = 'L': LOWER TRIANGLE OF A CONTAINS ELEMENTARY REFLECTORS * FROM DSYTRD. * * TRANS (INPUT) CHARACTER*1 * = 'N': NO TRANSPOSE, APPLY Q; * = 'T': TRANSPOSE, APPLY Q**T. * * M (INPUT) INTEGER * THE NUMBER OF ROWS OF THE MATRIX C. M >= 0. * * N (INPUT) INTEGER * THE NUMBER OF COLUMNS OF THE MATRIX C. N >= 0. * * A (INPUT) DOUBLE PRECISION ARRAY, DIMENSION * (LDA,M) IF SIDE = 'L' * (LDA,N) IF SIDE = 'R' * THE VECTORS WHICH DEFINE THE ELEMENTARY REFLECTORS, AS * RETURNED BY DSYTRD. * * LDA (INPUT) INTEGER * THE LEADING DIMENSION OF THE ARRAY A. * LDA >= MAX(1,M) IF SIDE = 'L'; LDA >= MAX(1,N) IF SIDE = 'R'. * * TAU (INPUT) DOUBLE PRECISION ARRAY, DIMENSION * (M-1) IF SIDE = 'L' * (N-1) IF SIDE = 'R' * TAU(I) MUST CONTAIN THE SCALAR FACTOR OF THE ELEMENTARY * REFLECTOR H(I), AS RETURNED BY DSYTRD. * * C (INPUT/OUTPUT) DOUBLE PRECISION ARRAY, DIMENSION (LDC,N) * ON ENTRY, THE M-BY-N MATRIX C. * ON EXIT, C IS OVERWRITTEN BY Q*C OR Q**T*C OR C*Q**T OR C*Q. * * LDC (INPUT) INTEGER * THE LEADING DIMENSION OF THE ARRAY C. LDC >= MAX(1,M). * * WORK (WORKSPACE) DOUBLE PRECISION ARRAY, DIMENSION (LWORK) * ON EXIT, IF INFO = 0, WORK(1) RETURNS THE OPTIMAL LWORK. * * LWORK (INPUT) INTEGER * THE DIMENSION OF THE ARRAY WORK. * IF SIDE = 'L', LWORK >= MAX(1,N); * IF SIDE = 'R', LWORK >= MAX(1,M). * FOR OPTIMUM PERFORMANCE LWORK >= N*NB IF SIDE = 'L', AND * LWORK >= M*NB IF SIDE = 'R', WHERE NB IS THE OPTIMAL * BLOCKSIZE. * * INFO (OUTPUT) INTEGER * = 0: SUCCESSFUL EXIT * < 0: IF INFO = -I, THE I-TH ARGUMENT HAD AN ILLEGAL VALUE * * ===================================================================== * * .. LOCAL SCALARS .. LOGICAL LEFT, UPPER INTEGER I1, I2, IINFO, MI, NI, NQ, NW * .. * .. EXTERNAL FUNCTIONS .. LOGICAL LSAME EXTERNAL LSAME * .. * .. EXTERNAL SUBROUTINES .. EXTERNAL DORMQL, DORMQR, XERBLA * .. * .. INTRINSIC FUNCTIONS .. INTRINSIC MAX * .. * .. EXECUTABLE STATEMENTS .. * * TEST THE INPUT ARGUMENTS * INFO = 0 LEFT = LSAME( SIDE, 'L' ) UPPER = LSAME( UPLO, 'U' ) * * NQ IS THE ORDER OF Q AND NW IS THE MINIMUM DIMENSION OF WORK * IF( LEFT ) THEN NQ = M NW = N ELSE NQ = N NW = M END IF IF( .NOT.LEFT .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN INFO = -1 ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN INFO = -2 ELSE IF( .NOT.LSAME( TRANS, 'N' ) .AND. .NOT.LSAME( TRANS, 'T' ) ) $ THEN INFO = -3 ELSE IF( M.LT.0 ) THEN INFO = -4 ELSE IF( N.LT.0 ) THEN INFO = -5 ELSE IF( LDA.LT.MAX( 1, NQ ) ) THEN INFO = -7 ELSE IF( LDC.LT.MAX( 1, M ) ) THEN INFO = -10 ELSE IF( LWORK.LT.MAX( 1, NW ) ) THEN INFO = -12 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DORMTR', -INFO ) RETURN END IF * * QUICK RETURN IF POSSIBLE * IF( M.EQ.0 .OR. N.EQ.0 .OR. NQ.EQ.1 ) THEN WORK( 1 ) = 1 RETURN END IF * IF( LEFT ) THEN MI = M - 1 NI = N ELSE MI = M NI = N - 1 END IF * IF( UPPER ) THEN * * Q WAS DETERMINED BY A CALL TO DSYTRD WITH UPLO = 'U' * CALL DORMQL( SIDE, TRANS, MI, NI, NQ-1, A( 1, 2 ), LDA, TAU, C, $ LDC, WORK, LWORK, IINFO ) ELSE * * Q WAS DETERMINED BY A CALL TO DSYTRD WITH UPLO = 'L' * IF( LEFT ) THEN I1 = 2 I2 = 1 ELSE I1 = 1 I2 = 2 END IF CALL DORMQR( SIDE, TRANS, MI, NI, NQ-1, A( 2, 1 ), LDA, TAU, $ C( I1, I2 ), LDC, WORK, LWORK, IINFO ) END IF RETURN * * END OF DORMTR * END CUT HERE............ CAT > DORMQR.F <<'CUT HERE............' SUBROUTINE DORMQR( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, $ WORK, LWORK, INFO ) * * -- LAPACK ROUTINE (VERSION 1.1) -- * UNIV. OF TENNESSEE, UNIV. OF CALIFORNIA BERKELEY, NAG LTD., * COURANT INSTITUTE, ARGONNE NATIONAL LAB, AND RICE UNIVERSITY * MARCH 31, 1993 * * .. SCALAR ARGUMENTS .. CHARACTER SIDE, TRANS INTEGER INFO, K, LDA, LDC, LWORK, M, N * .. * .. ARRAY ARGUMENTS .. DOUBLE PRECISION A( LDA, * ), C( LDC, * ), TAU( * ), $ WORK( LWORK ) * .. * * PURPOSE * ======= * * DORMQR OVERWRITES THE GENERAL REAL M-BY-N MATRIX C WITH * * SIDE = 'L' SIDE = 'R' * TRANS = 'N': Q * C C * Q * TRANS = 'T': Q**T * C C * Q**T * * WHERE Q IS A REAL ORTHOGONAL MATRIX DEFINED AS THE PRODUCT OF K * ELEMENTARY REFLECTORS * * Q = H(1) H(2) . . . H(K) * * AS RETURNED BY DGEQRF. Q IS OF ORDER M IF SIDE = 'L' AND OF ORDER N * IF SIDE = 'R'. * * ARGUMENTS * ========= * * SIDE (INPUT) CHARACTER*1 * = 'L': APPLY Q OR Q**T FROM THE LEFT; * = 'R': APPLY Q OR Q**T FROM THE RIGHT. * * TRANS (INPUT) CHARACTER*1 * = 'N': NO TRANSPOSE, APPLY Q; * = 'T': TRANSPOSE, APPLY Q**T. * * M (INPUT) INTEGER * THE NUMBER OF ROWS OF THE MATRIX C. M >= 0. * * N (INPUT) INTEGER * THE NUMBER OF COLUMNS OF THE MATRIX C. N >= 0. * * K (INPUT) INTEGER * THE NUMBER OF ELEMENTARY REFLECTORS WHOSE PRODUCT DEFINES * THE MATRIX Q. * IF SIDE = 'L', M >= K >= 0; * IF SIDE = 'R', N >= K >= 0. * * A (INPUT) DOUBLE PRECISION ARRAY, DIMENSION (LDA,K) * THE I-TH COLUMN MUST CONTAIN THE VECTOR WHICH DEFINES THE * ELEMENTARY REFLECTOR H(I), FOR I = 1,2,...,K, AS RETURNED BY * DGEQRF IN THE FIRST K COLUMNS OF ITS ARRAY ARGUMENT A. * A IS MODIFIED BY THE ROUTINE BUT RESTORED ON EXIT. * * LDA (INPUT) INTEGER * THE LEADING DIMENSION OF THE ARRAY A. * IF SIDE = 'L', LDA >= MAX(1,M); * IF SIDE = 'R', LDA >= MAX(1,N). * * TAU (INPUT) DOUBLE PRECISION ARRAY, DIMENSION (K) * TAU(I) MUST CONTAIN THE SCALAR FACTOR OF THE ELEMENTARY * REFLECTOR H(I), AS RETURNED BY DGEQRF. * * C (INPUT/OUTPUT) DOUBLE PRECISION ARRAY, DIMENSION (LDC,N) * ON ENTRY, THE M-BY-N MATRIX C. * ON EXIT, C IS OVERWRITTEN BY Q*C OR Q**T*C OR C*Q**T OR C*Q. * * LDC (INPUT) INTEGER * THE LEADING DIMENSION OF THE ARRAY C. LDC >= MAX(1,M). * * WORK (WORKSPACE) DOUBLE PRECISION ARRAY, DIMENSION (LWORK) * ON EXIT, IF INFO = 0, WORK(1) RETURNS THE OPTIMAL LWORK. * * LWORK (INPUT) INTEGER * THE DIMENSION OF THE ARRAY WORK. * IF SIDE = 'L', LWORK >= MAX(1,N); * IF SIDE = 'R', LWORK >= MAX(1,M). * FOR OPTIMUM PERFORMANCE LWORK >= N*NB IF SIDE = 'L', AND * LWORK >= M*NB IF SIDE = 'R', WHERE NB IS THE OPTIMAL * BLOCKSIZE. * * INFO (OUTPUT) INTEGER * = 0: SUCCESSFUL EXIT * < 0: IF INFO = -I, THE I-TH ARGUMENT HAD AN ILLEGAL VALUE * * ===================================================================== * * .. PARAMETERS .. INTEGER NBMAX, LDT PARAMETER ( NBMAX = 64, LDT = NBMAX+1 ) * .. * .. LOCAL SCALARS .. LOGICAL LEFT, NOTRAN INTEGER I, I1, I2, I3, IB, IC, IINFO, IWS, JC, LDWORK, $ MI, NB, NBMIN, NI, NQ, NW * .. * .. LOCAL ARRAYS .. DOUBLE PRECISION T( LDT, NBMAX ) * .. * .. EXTERNAL FUNCTIONS .. LOGICAL LSAME INTEGER ILAENV EXTERNAL LSAME, ILAENV * .. * .. EXTERNAL SUBROUTINES .. EXTERNAL DLARFB, DLARFT, DORM2R, XERBLA * .. * .. INTRINSIC FUNCTIONS .. INTRINSIC MAX, MIN * .. * .. EXECUTABLE STATEMENTS .. * * TEST THE INPUT ARGUMENTS * INFO = 0 LEFT = LSAME( SIDE, 'L' ) NOTRAN = LSAME( TRANS, 'N' ) * * NQ IS THE ORDER OF Q AND NW IS THE MINIMUM DIMENSION OF WORK * IF( LEFT ) THEN NQ = M NW = N ELSE NQ = N NW = M END IF IF( .NOT.LEFT .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN INFO = -1 ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN INFO = -2 ELSE IF( M.LT.0 ) THEN INFO = -3 ELSE IF( N.LT.0 ) THEN INFO = -4 ELSE IF( K.LT.0 .OR. K.GT.NQ ) THEN INFO = -5 ELSE IF( LDA.LT.MAX( 1, NQ ) ) THEN INFO = -7 ELSE IF( LDC.LT.MAX( 1, M ) ) THEN INFO = -10 ELSE IF( LWORK.LT.MAX( 1, NW ) ) THEN INFO = -12 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DORMQR', -INFO ) RETURN END IF * * QUICK RETURN IF POSSIBLE * IF( M.EQ.0 .OR. N.EQ.0 .OR. K.EQ.0 ) THEN WORK( 1 ) = 1 RETURN END IF * * DETERMINE THE BLOCK SIZE. NB MAY BE AT MOST NBMAX, WHERE NBMAX * IS USED TO DEFINE THE LOCAL ARRAY T. * NB = MIN( NBMAX, ILAENV( 1, 'DORMQR', SIDE // TRANS, M, N, K, $ -1 ) ) NBMIN = 2 LDWORK = NW IF( NB.GT.1 .AND. NB.LT.K ) THEN IWS = NW*NB IF( LWORK.LT.IWS ) THEN NB = LWORK / LDWORK NBMIN = MAX( 2, ILAENV( 2, 'DORMQR', SIDE // TRANS, M, N, K, $ -1 ) ) END IF ELSE IWS = NW END IF * IF( NB.LT.NBMIN .OR. NB.GE.K ) THEN * * USE UNBLOCKED CODE * CALL DORM2R( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, $ IINFO ) ELSE * * USE BLOCKED CODE * IF( ( LEFT .AND. .NOT.NOTRAN ) .OR. $ ( .NOT.LEFT .AND. NOTRAN ) ) THEN I1 = 1 I2 = K I3 = NB ELSE I1 = ( ( K-1 ) / NB )*NB + 1 I2 = 1 I3 = -NB END IF * IF( LEFT ) THEN NI = N JC = 1 ELSE MI = M IC = 1 END IF * DO 10 I = I1, I2, I3 IB = MIN( NB, K-I+1 ) * * FORM THE TRIANGULAR FACTOR OF THE BLOCK REFLECTOR * H = H(I) H(I+1) . . . H(I+IB-1) * CALL DLARFT( 'FORWARD', 'COLUMNWISE', NQ-I+1, IB, A( I, I ), $ LDA, TAU( I ), T, LDT ) IF( LEFT ) THEN * * H OR H' IS APPLIED TO C(I:M,1:N) * MI = M - I + 1 IC = I ELSE * * H OR H' IS APPLIED TO C(1:M,I:N) * NI = N - I + 1 JC = I END IF * * APPLY H OR H' * CALL DLARFB( SIDE, TRANS, 'FORWARD', 'COLUMNWISE', MI, NI, $ IB, A( I, I ), LDA, T, LDT, C( IC, JC ), LDC, $ WORK, LDWORK ) 10 CONTINUE END IF WORK( 1 ) = IWS RETURN * * END OF DORMQR * END CUT HERE............ CAT > DORM2R.F <<'CUT HERE............' SUBROUTINE DORM2R( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, $ WORK, INFO ) * * -- LAPACK ROUTINE (VERSION 1.1) -- * UNIV. OF TENNESSEE, UNIV. OF CALIFORNIA BERKELEY, NAG LTD., * COURANT INSTITUTE, ARGONNE NATIONAL LAB, AND RICE UNIVERSITY * FEBRUARY 29, 1992 * * .. SCALAR ARGUMENTS .. CHARACTER SIDE, TRANS INTEGER INFO, K, LDA, LDC, M, N * .. * .. ARRAY ARGUMENTS .. DOUBLE PRECISION A( LDA, * ), C( LDC, * ), TAU( * ), WORK( * ) * .. * * PURPOSE * ======= * * DORM2R OVERWRITES THE GENERAL REAL M BY N MATRIX C WITH * * Q * C IF SIDE = 'L' AND TRANS = 'N', OR * * Q'* C IF SIDE = 'L' AND TRANS = 'T', OR * * C * Q IF SIDE = 'R' AND TRANS = 'N', OR * * C * Q' IF SIDE = 'R' AND TRANS = 'T', * * WHERE Q IS A REAL ORTHOGONAL MATRIX DEFINED AS THE PRODUCT OF K * ELEMENTARY REFLECTORS * * Q = H(1) H(2) . . . H(K) * * AS RETURNED BY DGEQRF. Q IS OF ORDER M IF SIDE = 'L' AND OF ORDER N * IF SIDE = 'R'. * * ARGUMENTS * ========= * * SIDE (INPUT) CHARACTER*1 * = 'L': APPLY Q OR Q' FROM THE LEFT * = 'R': APPLY Q OR Q' FROM THE RIGHT * * TRANS (INPUT) CHARACTER*1 * = 'N': APPLY Q (NO TRANSPOSE) * = 'T': APPLY Q' (TRANSPOSE) * * M (INPUT) INTEGER * THE NUMBER OF ROWS OF THE MATRIX C. M >= 0. * * N (INPUT) INTEGER * THE NUMBER OF COLUMNS OF THE MATRIX C. N >= 0. * * K (INPUT) INTEGER * THE NUMBER OF ELEMENTARY REFLECTORS WHOSE PRODUCT DEFINES * THE MATRIX Q. * IF SIDE = 'L', M >= K >= 0; * IF SIDE = 'R', N >= K >= 0. * * A (INPUT) DOUBLE PRECISION ARRAY, DIMENSION (LDA,K) * THE I-TH COLUMN MUST CONTAIN THE VECTOR WHICH DEFINES THE * ELEMENTARY REFLECTOR H(I), FOR I = 1,2,...,K, AS RETURNED BY * DGEQRF IN THE FIRST K COLUMNS OF ITS ARRAY ARGUMENT A. * A IS MODIFIED BY THE ROUTINE BUT RESTORED ON EXIT. * * LDA (INPUT) INTEGER * THE LEADING DIMENSION OF THE ARRAY A. * IF SIDE = 'L', LDA >= MAX(1,M); * IF SIDE = 'R', LDA >= MAX(1,N). * * TAU (INPUT) DOUBLE PRECISION ARRAY, DIMENSION (K) * TAU(I) MUST CONTAIN THE SCALAR FACTOR OF THE ELEMENTARY * REFLECTOR H(I), AS RETURNED BY DGEQRF. * * C (INPUT/OUTPUT) DOUBLE PRECISION ARRAY, DIMENSION (LDC,N) * ON ENTRY, THE M BY N MATRIX C. * ON EXIT, C IS OVERWRITTEN BY Q*C OR Q'*C OR C*Q' OR C*Q. * * LDC (INPUT) INTEGER * THE LEADING DIMENSION OF THE ARRAY C. LDC >= MAX(1,M). * * WORK (WORKSPACE) DOUBLE PRECISION ARRAY, DIMENSION * (N) IF SIDE = 'L', * (M) IF SIDE = 'R' * * INFO (OUTPUT) INTEGER * = 0: SUCCESSFUL EXIT * < 0: IF INFO = -I, THE I-TH ARGUMENT HAD AN ILLEGAL VALUE * * ===================================================================== * * .. PARAMETERS .. DOUBLE PRECISION ONE PARAMETER ( ONE = 1.0D+0 ) * .. * .. LOCAL SCALARS .. LOGICAL LEFT, NOTRAN INTEGER I, I1, I2, I3, IC, JC, MI, NI, NQ DOUBLE PRECISION AII * .. * .. EXTERNAL FUNCTIONS .. LOGICAL LSAME EXTERNAL LSAME * .. * .. EXTERNAL SUBROUTINES .. EXTERNAL DLARF, XERBLA * .. * .. INTRINSIC FUNCTIONS .. INTRINSIC MAX * .. * .. EXECUTABLE STATEMENTS .. * * TEST THE INPUT ARGUMENTS * INFO = 0 LEFT = LSAME( SIDE, 'L' ) NOTRAN = LSAME( TRANS, 'N' ) * * NQ IS THE ORDER OF Q * IF( LEFT ) THEN NQ = M ELSE NQ = N END IF IF( .NOT.LEFT .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN INFO = -1 ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN INFO = -2 ELSE IF( M.LT.0 ) THEN INFO = -3 ELSE IF( N.LT.0 ) THEN INFO = -4 ELSE IF( K.LT.0 .OR. K.GT.NQ ) THEN INFO = -5 ELSE IF( LDA.LT.MAX( 1, NQ ) ) THEN INFO = -7 ELSE IF( LDC.LT.MAX( 1, M ) ) THEN INFO = -10 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DORM2R', -INFO ) RETURN END IF * * QUICK RETURN IF POSSIBLE * IF( M.EQ.0 .OR. N.EQ.0 .OR. K.EQ.0 ) $ RETURN * IF( ( LEFT .AND. .NOT.NOTRAN ) .OR. ( .NOT.LEFT .AND. NOTRAN ) ) $ THEN I1 = 1 I2 = K I3 = 1 ELSE I1 = K I2 = 1 I3 = -1 END IF * IF( LEFT ) THEN NI = N JC = 1 ELSE MI = M IC = 1 END IF * DO 10 I = I1, I2, I3 IF( LEFT ) THEN * * H(I) IS APPLIED TO C(I:M,1:N) * MI = M - I + 1 IC = I ELSE * * H(I) IS APPLIED TO C(1:M,I:N) * NI = N - I + 1 JC = I END IF * * APPLY H(I) * AII = A( I, I ) A( I, I ) = ONE CALL DLARF( SIDE, MI, NI, A( I, I ), 1, TAU( I ), C( IC, JC ), $ LDC, WORK ) A( I, I ) = AII 10 CONTINUE RETURN * * END OF DORM2R * END CUT HERE............ CAT > DORMQL.F <<'CUT HERE............' SUBROUTINE DORMQL( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, $ WORK, LWORK, INFO ) * * -- LAPACK ROUTINE (VERSION 1.1) -- * UNIV. OF TENNESSEE, UNIV. OF CALIFORNIA BERKELEY, NAG LTD., * COURANT INSTITUTE, ARGONNE NATIONAL LAB, AND RICE UNIVERSITY * MARCH 31, 1993 * * .. SCALAR ARGUMENTS .. CHARACTER SIDE, TRANS INTEGER INFO, K, LDA, LDC, LWORK, M, N * .. * .. ARRAY ARGUMENTS .. DOUBLE PRECISION A( LDA, * ), C( LDC, * ), TAU( * ), $ WORK( LWORK ) * .. * * PURPOSE * ======= * * DORMQL OVERWRITES THE GENERAL REAL M-BY-N MATRIX C WITH * * SIDE = 'L' SIDE = 'R' * TRANS = 'N': Q * C C * Q * TRANS = 'T': Q**T * C C * Q**T * * WHERE Q IS A REAL ORTHOGONAL MATRIX DEFINED AS THE PRODUCT OF K * ELEMENTARY REFLECTORS * * Q = H(K) . . . H(2) H(1) * * AS RETURNED BY DGEQLF. Q IS OF ORDER M IF SIDE = 'L' AND OF ORDER N * IF SIDE = 'R'. * * ARGUMENTS * ========= * * SIDE (INPUT) CHARACTER*1 * = 'L': APPLY Q OR Q**T FROM THE LEFT; * = 'R': APPLY Q OR Q**T FROM THE RIGHT. * * TRANS (INPUT) CHARACTER*1 * = 'N': NO TRANSPOSE, APPLY Q; * = 'T': TRANSPOSE, APPLY Q**T. * * M (INPUT) INTEGER * THE NUMBER OF ROWS OF THE MATRIX C. M >= 0. * * N (INPUT) INTEGER * THE NUMBER OF COLUMNS OF THE MATRIX C. N >= 0. * * K (INPUT) INTEGER * THE NUMBER OF ELEMENTARY REFLECTORS WHOSE PRODUCT DEFINES * THE MATRIX Q. * IF SIDE = 'L', M >= K >= 0; * IF SIDE = 'R', N >= K >= 0. * * A (INPUT) DOUBLE PRECISION ARRAY, DIMENSION (LDA,K) * THE I-TH COLUMN MUST CONTAIN THE VECTOR WHICH DEFINES THE * ELEMENTARY REFLECTOR H(I), FOR I = 1,2,...,K, AS RETURNED BY * DGEQLF IN THE LAST K COLUMNS OF ITS ARRAY ARGUMENT A. * A IS MODIFIED BY THE ROUTINE BUT RESTORED ON EXIT. * * LDA (INPUT) INTEGER * THE LEADING DIMENSION OF THE ARRAY A. * IF SIDE = 'L', LDA >= MAX(1,M); * IF SIDE = 'R', LDA >= MAX(1,N). * * TAU (INPUT) DOUBLE PRECISION ARRAY, DIMENSION (K) * TAU(I) MUST CONTAIN THE SCALAR FACTOR OF THE ELEMENTARY * REFLECTOR H(I), AS RETURNED BY DGEQLF. * * C (INPUT/OUTPUT) DOUBLE PRECISION ARRAY, DIMENSION (LDC,N) * ON ENTRY, THE M-BY-N MATRIX C. * ON EXIT, C IS OVERWRITTEN BY Q*C OR Q**T*C OR C*Q**T OR C*Q. * * LDC (INPUT) INTEGER * THE LEADING DIMENSION OF THE ARRAY C. LDC >= MAX(1,M). * * WORK (WORKSPACE) DOUBLE PRECISION ARRAY, DIMENSION (LWORK) * ON EXIT, IF INFO = 0, WORK(1) RETURNS THE OPTIMAL LWORK. * * LWORK (INPUT) INTEGER * THE DIMENSION OF THE ARRAY WORK. * IF SIDE = 'L', LWORK >= MAX(1,N); * IF SIDE = 'R', LWORK >= MAX(1,M). * FOR OPTIMUM PERFORMANCE LWORK >= N*NB IF SIDE = 'L', AND * LWORK >= M*NB IF SIDE = 'R', WHERE NB IS THE OPTIMAL * BLOCKSIZE. * * INFO (OUTPUT) INTEGER * = 0: SUCCESSFUL EXIT * < 0: IF INFO = -I, THE I-TH ARGUMENT HAD AN ILLEGAL VALUE * * ===================================================================== * * .. PARAMETERS .. INTEGER NBMAX, LDT PARAMETER ( NBMAX = 64, LDT = NBMAX+1 ) * .. * .. LOCAL SCALARS .. LOGICAL LEFT, NOTRAN INTEGER I, I1, I2, I3, IB, IINFO, IWS, LDWORK, MI, NB, $ NBMIN, NI, NQ, NW * .. * .. LOCAL ARRAYS .. DOUBLE PRECISION T( LDT, NBMAX ) * .. * .. EXTERNAL FUNCTIONS .. LOGICAL LSAME INTEGER ILAENV EXTERNAL LSAME, ILAENV * .. * .. EXTERNAL SUBROUTINES .. EXTERNAL DLARFB, DLARFT, DORM2L, XERBLA * .. * .. INTRINSIC FUNCTIONS .. INTRINSIC MAX, MIN * .. * .. EXECUTABLE STATEMENTS .. * * TEST THE INPUT ARGUMENTS * INFO = 0 LEFT = LSAME( SIDE, 'L' ) NOTRAN = LSAME( TRANS, 'N' ) * * NQ IS THE ORDER OF Q AND NW IS THE MINIMUM DIMENSION OF WORK * IF( LEFT ) THEN NQ = M NW = N ELSE NQ = N NW = M END IF IF( .NOT.LEFT .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN INFO = -1 ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN INFO = -2 ELSE IF( M.LT.0 ) THEN INFO = -3 ELSE IF( N.LT.0 ) THEN INFO = -4 ELSE IF( K.LT.0 .OR. K.GT.NQ ) THEN INFO = -5 ELSE IF( LDA.LT.MAX( 1, NQ ) ) THEN INFO = -7 ELSE IF( LDC.LT.MAX( 1, M ) ) THEN INFO = -10 ELSE IF( LWORK.LT.MAX( 1, NW ) ) THEN INFO = -12 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DORMQL', -INFO ) RETURN END IF * * QUICK RETURN IF POSSIBLE * IF( M.EQ.0 .OR. N.EQ.0 .OR. K.EQ.0 ) THEN WORK( 1 ) = 1 RETURN END IF * * DETERMINE THE BLOCK SIZE. NB MAY BE AT MOST NBMAX, WHERE NBMAX * IS USED TO DEFINE THE LOCAL ARRAY T. * NB = MIN( NBMAX, ILAENV( 1, 'DORMQL', SIDE // TRANS, M, N, K, $ -1 ) ) NBMIN = 2 LDWORK = NW IF( NB.GT.1 .AND. NB.LT.K ) THEN IWS = NW*NB IF( LWORK.LT.IWS ) THEN NB = LWORK / LDWORK NBMIN = MAX( 2, ILAENV( 2, 'DORMQL', SIDE // TRANS, M, N, K, $ -1 ) ) END IF ELSE IWS = NW END IF * IF( NB.LT.NBMIN .OR. NB.GE.K ) THEN * * USE UNBLOCKED CODE * CALL DORM2L( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, $ IINFO ) ELSE * * USE BLOCKED CODE * IF( ( LEFT .AND. NOTRAN ) .OR. $ ( .NOT.LEFT .AND. .NOT.NOTRAN ) ) THEN I1 = 1 I2 = K I3 = NB ELSE I1 = ( ( K-1 ) / NB )*NB + 1 I2 = 1 I3 = -NB END IF * IF( LEFT ) THEN NI = N ELSE MI = M END IF * DO 10 I = I1, I2, I3 IB = MIN( NB, K-I+1 ) * * FORM THE TRIANGULAR FACTOR OF THE BLOCK REFLECTOR * H = H(I+IB-1) . . . H(I+1) H(I) * CALL DLARFT( 'BACKWARD', 'COLUMNWISE', NQ-K+I+IB-1, IB, $ A( 1, I ), LDA, TAU( I ), T, LDT ) IF( LEFT ) THEN * * H OR H' IS APPLIED TO C(1:M-K+I+IB-1,1:N) * MI = M - K + I + IB - 1 ELSE * * H OR H' IS APPLIED TO C(1:M,1:N-K+I+IB-1) * NI = N - K + I + IB - 1 END IF * * APPLY H OR H' * CALL DLARFB( SIDE, TRANS, 'BACKWARD', 'COLUMNWISE', MI, NI, $ IB, A( 1, I ), LDA, T, LDT, C, LDC, WORK, $ LDWORK ) 10 CONTINUE END IF WORK( 1 ) = IWS RETURN * * END OF DORMQL * END CUT HERE............ CAT > DORM2L.F <<'CUT HERE............' SUBROUTINE DORM2L( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, $ WORK, INFO ) * * -- LAPACK ROUTINE (VERSION 1.1) -- * UNIV. OF TENNESSEE, UNIV. OF CALIFORNIA BERKELEY, NAG LTD., * COURANT INSTITUTE, ARGONNE NATIONAL LAB, AND RICE UNIVERSITY * FEBRUARY 29, 1992 * * .. SCALAR ARGUMENTS .. CHARACTER SIDE, TRANS INTEGER INFO, K, LDA, LDC, M, N * .. * .. ARRAY ARGUMENTS .. DOUBLE PRECISION A( LDA, * ), C( LDC, * ), TAU( * ), WORK( * ) * .. * * PURPOSE * ======= * * DORM2L OVERWRITES THE GENERAL REAL M BY N MATRIX C WITH * * Q * C IF SIDE = 'L' AND TRANS = 'N', OR * * Q'* C IF SIDE = 'L' AND TRANS = 'T', OR * * C * Q IF SIDE = 'R' AND TRANS = 'N', OR * * C * Q' IF SIDE = 'R' AND TRANS = 'T', * * WHERE Q IS A REAL ORTHOGONAL MATRIX DEFINED AS THE PRODUCT OF K * ELEMENTARY REFLECTORS * * Q = H(K) . . . H(2) H(1) * * AS RETURNED BY DGEQLF. Q IS OF ORDER M IF SIDE = 'L' AND OF ORDER N * IF SIDE = 'R'. * * ARGUMENTS * ========= * * SIDE (INPUT) CHARACTER*1 * = 'L': APPLY Q OR Q' FROM THE LEFT * = 'R': APPLY Q OR Q' FROM THE RIGHT * * TRANS (INPUT) CHARACTER*1 * = 'N': APPLY Q (NO TRANSPOSE) * = 'T': APPLY Q' (TRANSPOSE) * * M (INPUT) INTEGER * THE NUMBER OF ROWS OF THE MATRIX C. M >= 0. * * N (INPUT) INTEGER * THE NUMBER OF COLUMNS OF THE MATRIX C. N >= 0. * * K (INPUT) INTEGER * THE NUMBER OF ELEMENTARY REFLECTORS WHOSE PRODUCT DEFINES * THE MATRIX Q. * IF SIDE = 'L', M >= K >= 0; * IF SIDE = 'R', N >= K >= 0. * * A (INPUT) DOUBLE PRECISION ARRAY, DIMENSION (LDA,K) * THE I-TH COLUMN MUST CONTAIN THE VECTOR WHICH DEFINES THE * ELEMENTARY REFLECTOR H(I), FOR I = 1,2,...,K, AS RETURNED BY * DGEQLF IN THE LAST K COLUMNS OF ITS ARRAY ARGUMENT A. * A IS MODIFIED BY THE ROUTINE BUT RESTORED ON EXIT. * * LDA (INPUT) INTEGER * THE LEADING DIMENSION OF THE ARRAY A. * IF SIDE = 'L', LDA >= MAX(1,M); * IF SIDE = 'R', LDA >= MAX(1,N). * * TAU (INPUT) DOUBLE PRECISION ARRAY, DIMENSION (K) * TAU(I) MUST CONTAIN THE SCALAR FACTOR OF THE ELEMENTARY * REFLECTOR H(I), AS RETURNED BY DGEQLF. * * C (INPUT/OUTPUT) DOUBLE PRECISION ARRAY, DIMENSION (LDC,N) * ON ENTRY, THE M BY N MATRIX C. * ON EXIT, C IS OVERWRITTEN BY Q*C OR Q'*C OR C*Q' OR C*Q. * * LDC (INPUT) INTEGER * THE LEADING DIMENSION OF THE ARRAY C. LDC >= MAX(1,M). * * WORK (WORKSPACE) DOUBLE PRECISION ARRAY, DIMENSION * (N) IF SIDE = 'L', * (M) IF SIDE = 'R' * * INFO (OUTPUT) INTEGER * = 0: SUCCESSFUL EXIT * < 0: IF INFO = -I, THE I-TH ARGUMENT HAD AN ILLEGAL VALUE * * ===================================================================== * * .. PARAMETERS .. DOUBLE PRECISION ONE PARAMETER ( ONE = 1.0D+0 ) * .. * .. LOCAL SCALARS .. LOGICAL LEFT, NOTRAN INTEGER I, I1, I2, I3, MI, NI, NQ DOUBLE PRECISION AII * .. * .. EXTERNAL FUNCTIONS .. LOGICAL LSAME EXTERNAL LSAME * .. * .. EXTERNAL SUBROUTINES .. EXTERNAL DLARF, XERBLA * .. * .. INTRINSIC FUNCTIONS .. INTRINSIC MAX * .. * .. EXECUTABLE STATEMENTS .. * * TEST THE INPUT ARGUMENTS * INFO = 0 LEFT = LSAME( SIDE, 'L' ) NOTRAN = LSAME( TRANS, 'N' ) * * NQ IS THE ORDER OF Q * IF( LEFT ) THEN NQ = M ELSE NQ = N END IF IF( .NOT.LEFT .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN INFO = -1 ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN INFO = -2 ELSE IF( M.LT.0 ) THEN INFO = -3 ELSE IF( N.LT.0 ) THEN INFO = -4 ELSE IF( K.LT.0 .OR. K.GT.NQ ) THEN INFO = -5 ELSE IF( LDA.LT.MAX( 1, NQ ) ) THEN INFO = -7 ELSE IF( LDC.LT.MAX( 1, M ) ) THEN INFO = -10 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DORM2L', -INFO ) RETURN END IF * * QUICK RETURN IF POSSIBLE * IF( M.EQ.0 .OR. N.EQ.0 .OR. K.EQ.0 ) $ RETURN * IF( ( LEFT .AND. NOTRAN ) .OR. ( .NOT.LEFT .AND. .NOT.NOTRAN ) ) $ THEN I1 = 1 I2 = K I3 = 1 ELSE I1 = K I2 = 1 I3 = -1 END IF * IF( LEFT ) THEN NI = N ELSE MI = M END IF * DO 10 I = I1, I2, I3 IF( LEFT ) THEN * * H(I) IS APPLIED TO C(1:M-K+I,1:N) * MI = M - K + I ELSE * * H(I) IS APPLIED TO C(1:M,1:N-K+I) * NI = N - K + I END IF * * APPLY H(I) * AII = A( NQ-K+I, I ) A( NQ-K+I, I ) = ONE CALL DLARF( SIDE, MI, NI, A( 1, I ), 1, TAU( I ), C, LDC, $ WORK ) A( NQ-K+I, I ) = AII 10 CONTINUE RETURN * * END OF DORM2L * END CUT HERE............ CAT > DSTEIN.F <<'CUT HERE............' SUBROUTINE DSTEIN( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, $ IWORK, IFAIL, INFO ) * * -- LAPACK ROUTINE (VERSION 1.1) -- * UNIV. OF TENNESSEE, UNIV. OF CALIFORNIA BERKELEY, NAG LTD., * COURANT INSTITUTE, ARGONNE NATIONAL LAB, AND RICE UNIVERSITY * MARCH 31, 1993 * * .. SCALAR ARGUMENTS .. INTEGER INFO, LDZ, M, N * .. * .. ARRAY ARGUMENTS .. INTEGER IBLOCK( * ), IFAIL( * ), ISPLIT( * ), $ IWORK( * ) DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * ) * .. * * PURPOSE * ======= * * DSTEIN COMPUTES THE EIGENVECTORS OF A REAL SYMMETRIC TRIDIAGONAL * MATRIX T CORRESPONDING TO SPECIFIED EIGENVALUES, USING INVERSE * ITERATION. * * THE MAXIMUM NUMBER OF ITERATIONS ALLOWED FOR EACH EIGENVECTOR IS * SPECIFIED BY AN INTERNAL PARAMETER MAXITS (CURRENTLY SET TO 5). * * ARGUMENTS * ========= * * N (INPUT) INTEGER * THE ORDER OF THE MATRIX. N >= 0. * * D (INPUT) DOUBLE PRECISION ARRAY, DIMENSION (N) * THE N DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX T. * * E (INPUT) DOUBLE PRECISION ARRAY, DIMENSION (N) * THE (N-1) SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX * T, IN ELEMENTS 1 TO N-1. E(N) NEED NOT BE SET. * * M (INPUT) INTEGER * THE NUMBER OF EIGENVECTORS TO BE FOUND. 0 <= M <= N. * * W (INPUT) DOUBLE PRECISION ARRAY, DIMENSION (N) * THE FIRST M ELEMENTS OF W CONTAIN THE EIGENVALUES FOR * WHICH EIGENVECTORS ARE TO BE COMPUTED. THE EIGENVALUES * SHOULD BE GROUPED BY SPLIT-OFF BLOCK AND ORDERED FROM * SMALLEST TO LARGEST WITHIN THE BLOCK. ( THE OUTPUT ARRAY * W FROM DSTEBZ WITH ORDER = 'B' IS EXPECTED HERE. ) * * IBLOCK (INPUT) INTEGER ARRAY, DIMENSION (N) * THE SUBMATRIX INDICES ASSOCIATED WITH THE CORRESPONDING * EIGENVALUES IN W; IBLOCK(I)=1 IF EIGENVALUE W(I) BELONGS TO * THE FIRST SUBMATRIX FROM THE TOP, =2 IF W(I) BELONGS TO * THE SECOND SUBMATRIX, ETC. ( THE OUTPUT ARRAY IBLOCK * FROM DSTEBZ IS EXPECTED HERE. ) * * ISPLIT (INPUT) INTEGER ARRAY, DIMENSION (N) * THE SPLITTING POINTS, AT WHICH T BREAKS UP INTO SUBMATRICES. * THE FIRST SUBMATRIX CONSISTS OF ROWS/COLUMNS 1 TO * ISPLIT( 1 ), THE SECOND OF ROWS/COLUMNS ISPLIT( 1 )+1 * THROUGH ISPLIT( 2 ), ETC. * ( THE OUTPUT ARRAY ISPLIT FROM DSTEBZ IS EXPECTED HERE. ) * * Z (OUTPUT) DOUBLE PRECISION ARRAY, DIMENSION (LDZ, M) * THE COMPUTED EIGENVECTORS. THE EIGENVECTOR ASSOCIATED * WITH THE EIGENVALUE W(I) IS STORED IN THE I-TH COLUMN OF * Z. ANY VECTOR WHICH FAILS TO CONVERGE IS SET TO ITS CURRENT * ITERATE AFTER MAXITS ITERATIONS. * * LDZ (INPUT) INTEGER * THE LEADING DIMENSION OF THE ARRAY Z. LDZ >= MAX(1,N). * * WORK (WORKSPACE) DOUBLE PRECISION ARRAY, DIMENSION (5*N) * * IWORK (WORKSPACE) INTEGER ARRAY, DIMENSION (N) * * IFAIL (OUTPUT) INTEGER ARRAY, DIMENSION (M) * ON NORMAL EXIT, ALL ELEMENTS OF IFAIL ARE ZERO. * IF ONE OR MORE EIGENVECTORS FAIL TO CONVERGE AFTER * MAXITS ITERATIONS, THEN THEIR INDICES ARE STORED IN * ARRAY IFAIL. * * INFO (OUTPUT) INTEGER * = 0: SUCCESSFUL EXIT. * < 0: IF INFO = -I, THE I-TH ARGUMENT HAD AN ILLEGAL VALUE * > 0: IF INFO = I, THEN I EIGENVECTORS FAILED TO CONVERGE * IN MAXITS ITERATIONS. THEIR INDICES ARE STORED IN * ARRAY IFAIL. * * INTERNAL PARAMETERS * =================== * * MAXITS INTEGER, DEFAULT = 5 * THE MAXIMUM NUMBER OF ITERATIONS PERFORMED. * * EXTRA INTEGER, DEFAULT = 2 * THE NUMBER OF ITERATIONS PERFORMED AFTER NORM GROWTH * CRITERION IS SATISFIED, SHOULD BE AT LEAST 1. * * ===================================================================== * * .. PARAMETERS .. DOUBLE PRECISION ZERO, ONE, TEN, ODM3, ODM1 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TEN = 1.0D1, $ ODM3 = 1.0D-3, ODM1 = 1.0D-1 ) INTEGER MAXITS, EXTRA PARAMETER ( MAXITS = 5, EXTRA = 2 ) * .. * .. LOCAL SCALARS .. INTEGER B1, BLKSIZ, BN, GPIND, I, IINFO, INDRV1, $ INDRV2, INDRV3, INDRV4, INDRV5, ITS, J, J1, $ JBLK, JMAX, NBLK, NRMCHK DOUBLE PRECISION DTPCRT, EPS, EPS1, NRM, ONENRM, ORTOL, PERTOL, $ SCL, SEP, TOL, XJ, XJM, ZTR * .. * .. LOCAL ARRAYS .. INTEGER ISEED( 4 ) * .. * .. EXTERNAL FUNCTIONS .. INTEGER IDAMAX DOUBLE PRECISION DASUM, DDOT, DLAMCH, DNRM2 EXTERNAL IDAMAX, DASUM, DDOT, DLAMCH, DNRM2 * .. * .. EXTERNAL SUBROUTINES .. EXTERNAL DAXPY, DCOPY, DLAGTF, DLAGTS, DLARNV, DSCAL, $ XERBLA * .. * .. INTRINSIC FUNCTIONS .. INTRINSIC ABS, MAX, SQRT * .. * .. EXECUTABLE STATEMENTS .. * * TEST THE INPUT PARAMETERS. * INFO = 0 DO 10 I = 1, M IFAIL( I ) = 0 10 CONTINUE * IF( N.LT.0 ) THEN INFO = -1 ELSE IF( M.LT.0 .OR. M.GT.N ) THEN INFO = -4 ELSE IF( LDZ.LT.MAX( 1, N ) ) THEN INFO = -9 ELSE DO 20 J = 2, M IF( IBLOCK( J ).LT.IBLOCK( J-1 ) ) THEN INFO = -6 GO TO 30 END IF IF( IBLOCK( J ).EQ.IBLOCK( J-1 ) .AND. W( J ).LT.W( J-1 ) ) $ THEN INFO = -5 GO TO 30 END IF 20 CONTINUE 30 CONTINUE END IF * IF( INFO.NE.0 ) THEN CALL XERBLA( 'DSTEIN', -INFO ) RETURN END IF * * QUICK RETURN IF POSSIBLE * IF( N.EQ.0 .OR. M.EQ.0 ) THEN RETURN ELSE IF( N.EQ.1 ) THEN Z( 1, 1 ) = ONE RETURN END IF * * GET MACHINE CONSTANTS. * EPS = DLAMCH( 'PRECISION' ) * * INITIALIZE SEED FOR RANDOM NUMBER GENERATOR DLARNV. * DO 40 I = 1, 4 ISEED( I ) = 1 40 CONTINUE * * INITIALIZE POINTERS. * INDRV1 = 0 INDRV2 = INDRV1 + N INDRV3 = INDRV2 + N INDRV4 = INDRV3 + N INDRV5 = INDRV4 + N * * COMPUTE EIGENVECTORS OF MATRIX BLOCKS. * GPIND = 1 J1 = 1 DO 160 NBLK = 1, IBLOCK( M ) * * FIND STARTING AND ENDING INDICES OF BLOCK NBLK. * IF( NBLK.EQ.1 ) THEN B1 = 1 ELSE B1 = ISPLIT( NBLK-1 ) + 1 END IF BN = ISPLIT( NBLK ) BLKSIZ = BN - B1 + 1 IF( BLKSIZ.EQ.1 ) $ GO TO 60 * * COMPUTE REORTHOGONALIZATION CRITERION AND STOPPING CRITERION. * ONENRM = ABS( D( B1 ) ) + ABS( E( B1 ) ) ONENRM = MAX( ONENRM, ABS( D( BN ) )+ABS( E( BN-1 ) ) ) DO 50 I = B1 + 1, BN - 1 ONENRM = MAX( ONENRM, ABS( D( I ) )+ABS( E( I-1 ) )+ $ ABS( E( I ) ) ) 50 CONTINUE ORTOL = ODM3*ONENRM * DTPCRT = SQRT( ODM1 / BLKSIZ ) * * LOOP THROUGH EIGENVALUES OF BLOCK NBLK. * 60 CONTINUE JBLK = 0 DO 150 J = J1, M IF( IBLOCK( J ).NE.NBLK ) THEN J1 = J GO TO 160 END IF JBLK = JBLK + 1 XJ = W( J ) * * SKIP ALL THE WORK IF THE BLOCK SIZE IS ONE. * IF( BLKSIZ.EQ.1 ) THEN WORK( INDRV1+1 ) = ONE GO TO 120 END IF * * IF EIGENVALUES J AND J-1 ARE TOO CLOSE, ADD A RELATIVELY * SMALL PERTURBATION. * IF( JBLK.GT.1 ) THEN EPS1 = ABS( EPS*XJ ) PERTOL = TEN*EPS1 SEP = XJ - XJM IF( SEP.LT.PERTOL ) $ XJ = XJM + PERTOL END IF * ITS = 0 NRMCHK = 0 * * GET RANDOM STARTING VECTOR. * CALL DLARNV( 3, ISEED, BLKSIZ, WORK( INDRV1+1 ) ) * * COPY THE MATRIX T SO IT WON'T BE DESTROYED IN FACTORIZATION. * CALL DCOPY( BLKSIZ, D( B1 ), 1, WORK( INDRV4+1 ), 1 ) CALL DCOPY( BLKSIZ-1, E( B1 ), 1, WORK( INDRV2+2 ), 1 ) CALL DCOPY( BLKSIZ-1, E( B1 ), 1, WORK( INDRV3+1 ), 1 ) * * COMPUTE LU FACTORS WITH PARTIAL PIVOTING ( PT = LU ) * TOL = ZERO CALL DLAGTF( BLKSIZ, WORK( INDRV4+1 ), XJ, WORK( INDRV2+2 ), $ WORK( INDRV3+1 ), TOL, WORK( INDRV5+1 ), IWORK, $ IINFO ) * * UPDATE ITERATION COUNT. * 70 CONTINUE ITS = ITS + 1 IF( ITS.GT.MAXITS ) $ GO TO 100 * * NORMALIZE AND SCALE THE RIGHTHAND SIDE VECTOR PB. * SCL = BLKSIZ*ONENRM*MAX( EPS, $ ABS( WORK( INDRV4+BLKSIZ ) ) ) / $ DASUM( BLKSIZ, WORK( INDRV1+1 ), 1 ) CALL DSCAL( BLKSIZ, SCL, WORK( INDRV1+1 ), 1 ) * * SOLVE THE SYSTEM LU = PB. * CALL DLAGTS( -1, BLKSIZ, WORK( INDRV4+1 ), WORK( INDRV2+2 ), $ WORK( INDRV3+1 ), WORK( INDRV5+1 ), IWORK, $ WORK( INDRV1+1 ), TOL, IINFO ) * * REORTHOGONALIZE BY MODIFIED GRAM-SCHMIDT IF EIGENVALUES ARE * CLOSE ENOUGH. * IF( JBLK.EQ.1 ) $ GO TO 90 IF( ABS( XJ-XJM ).GT.ORTOL ) $ GPIND = J IF( GPIND.NE.J ) THEN DO 80 I = GPIND, J - 1 ZTR = -DDOT( BLKSIZ, WORK( INDRV1+1 ), 1, Z( B1, I ), $ 1 ) CALL DAXPY( BLKSIZ, ZTR, Z( B1, I ), 1, $ WORK( INDRV1+1 ), 1 ) 80 CONTINUE END IF * * CHECK THE INFINITY NORM OF THE ITERATE. * 90 CONTINUE JMAX = IDAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 ) NRM = ABS( WORK( INDRV1+JMAX ) ) * * CONTINUE FOR ADDITIONAL ITERATIONS AFTER NORM REACHES * STOPPING CRITERION. * IF( NRM.LT.DTPCRT ) $ GO TO 70 NRMCHK = NRMCHK + 1 IF( NRMCHK.LT.EXTRA+1 ) $ GO TO 70 * GO TO 110 * * IF STOPPING CRITERION WAS NOT SATISFIED, UPDATE INFO AND * STORE EIGENVECTOR NUMBER IN ARRAY IFAIL. * 100 CONTINUE INFO = INFO + 1 IFAIL( INFO ) = J * * ACCEPT ITERATE AS JTH EIGENVECTOR. * 110 CONTINUE SCL = ONE / DNRM2( BLKSIZ, WORK( INDRV1+1 ), 1 ) JMAX = IDAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 ) IF( WORK( INDRV1+JMAX ).LT.ZERO ) $ SCL = -SCL CALL DSCAL( BLKSIZ, SCL, WORK( INDRV1+1 ), 1 ) 120 CONTINUE DO 130 I = 1, N Z( I, J ) = ZERO 130 CONTINUE DO 140 I = 1, BLKSIZ Z( B1+I-1, J ) = WORK( INDRV1+I ) 140 CONTINUE * * SAVE THE SHIFT TO CHECK EIGENVALUE SPACING AT NEXT * ITERATION. * XJM = XJ * 150 CONTINUE 160 CONTINUE * RETURN * * END OF DSTEIN * END CUT HERE............ CAT > DLAGTS.F <<'CUT HERE............' SUBROUTINE DLAGTS( JOB, N, A, B, C, D, IN, Y, TOL, INFO ) * * -- LAPACK AUXILIARY ROUTINE (VERSION 1.1) -- * UNIV. OF TENNESSEE, UNIV. OF CALIFORNIA BERKELEY, NAG LTD., * COURANT INSTITUTE, ARGONNE NATIONAL LAB, AND RICE UNIVERSITY * OCTOBER 31, 1992 * * .. SCALAR ARGUMENTS .. INTEGER INFO, JOB, N DOUBLE PRECISION TOL * .. * .. ARRAY ARGUMENTS .. INTEGER IN( * ) DOUBLE PRECISION A( * ), B( * ), C( * ), D( * ), Y( * ) * .. * * PURPOSE * ======= * * DLAGTS MAY BE USED TO SOLVE ONE OF THE SYSTEMS OF EQUATIONS * * (T - LAMBDA*I)*X = Y OR (T - LAMBDA*I)'*X = Y, * * WHERE T IS AN N BY N TRIDIAGONAL MATRIX, FOR X, FOLLOWING THE * FACTORIZATION OF (T - LAMBDA*I) AS * * (T - LAMBDA*I) = P*L*U , * * BY ROUTINE DLAGTF. THE CHOICE OF EQUATION TO BE SOLVED IS * CONTROLLED BY THE ARGUMENT JOB, AND IN EACH CASE THERE IS AN OPTION * TO PERTURB ZERO OR VERY SMALL DIAGONAL ELEMENTS OF U, THIS OPTION * BEING INTENDED FOR USE IN APPLICATIONS SUCH AS INVERSE ITERATION. * * ARGUMENTS * ========= * * JOB (INPUT) INTEGER * SPECIFIES THE JOB TO BE PERFORMED BY DLAGTS AS FOLLOWS: * = 1: THE EQUATIONS (T - LAMBDA*I)X = Y ARE TO BE SOLVED, * BUT DIAGONAL ELEMENTS OF U ARE NOT TO BE PERTURBED. * = -1: THE EQUATIONS (T - LAMBDA*I)X = Y ARE TO BE SOLVED * AND, IF OVERFLOW WOULD OTHERWISE OCCUR, THE DIAGONAL * ELEMENTS OF U ARE TO BE PERTURBED. SEE ARGUMENT TOL * BELOW. * = 2: THE EQUATIONS (T - LAMBDA*I)'X = Y ARE TO BE SOLVED, * BUT DIAGONAL ELEMENTS OF U ARE NOT TO BE PERTURBED. * = -2: THE EQUATIONS (T - LAMBDA*I)'X = Y ARE TO BE SOLVED * AND, IF OVERFLOW WOULD OTHERWISE OCCUR, THE DIAGONAL * ELEMENTS OF U ARE TO BE PERTURBED. SEE ARGUMENT TOL * BELOW. * * N (INPUT) INTEGER * THE ORDER OF THE MATRIX T. * * A (INPUT) DOUBLE PRECISION ARRAY, DIMENSION (N) * ON ENTRY, A MUST CONTAIN THE DIAGONAL ELEMENTS OF U AS * RETURNED FROM DLAGTF. * * B (INPUT) DOUBLE PRECISION ARRAY, DIMENSION (N-1) * ON ENTRY, B MUST CONTAIN THE FIRST SUPER-DIAGONAL ELEMENTS OF * U AS RETURNED FROM DLAGTF. * * C (INPUT) DOUBLE PRECISION ARRAY, DIMENSION (N-1) * ON ENTRY, C MUST CONTAIN THE SUB-DIAGONAL ELEMENTS OF L AS * RETURNED FROM DLAGTF. * * D (INPUT) DOUBLE PRECISION ARRAY, DIMENSION (N-2) * ON ENTRY, D MUST CONTAIN THE SECOND SUPER-DIAGONAL ELEMENTS * OF U AS RETURNED FROM DLAGTF. * * IN (INPUT) INTEGER ARRAY, DIMENSION (N) * ON ENTRY, IN MUST CONTAIN DETAILS OF THE MATRIX P AS RETURNED * FROM DLAGTF. * * Y (INPUT/OUTPUT) DOUBLE PRECISION ARRAY, DIMENSION (N) * ON ENTRY, THE RIGHT HAND SIDE VECTOR Y. * * ON EXIT, Y IS OVERWRITTEN BY THE SOLUTION VECTOR X. * * TOL (INPUT/OUTPUT) DOUBLE PRECISION * ON ENTRY WITH JOB .LT. 0, TOL SHOULD BE THE MINIMUM * PERTURBATION TO BE MADE TO VERY SMALL DIAGONAL ELEMENTS OF U. * TOL SHOULD NORMALLY BE CHOSEN AS ABOUT EPS*NORM(U), WHERE EPS * IS THE RELATIVE MACHINE PRECISION, BUT IF TOL IS SUPPLIED AS * NON-POSITIVE, THEN IT IS RESET TO EPS*MAX( ABS( U(I,J) ) ). * IF JOB .GT. 0 THEN TOL IS NOT REFERENCED. * * ON EXIT, TOL IS CHANGED AS DESCRIBED ABOVE, ONLY IF TOL IS * NON-POSITIVE ON ENTRY. OTHERWISE TOL IS UNCHANGED. * * INFO (OUTPUT) * = 0 : SUCCESSFUL EXIT * .LT. 0: IF INFO = -K, THE KTH ARGUMENT HAD AN ILLEGAL VALUE * .GT. 0: OVERFLOW WOULD OCCUR WHEN COMPUTING THE INFO(TH) * ELEMENT OF THE SOLUTION VECTOR X. THIS CAN ONLY OCCUR * WHEN JOB IS SUPPLIED AS POSITIVE AND EITHER MEANS * THAT A DIAGONAL ELEMENT OF U IS VERY SMALL, OR THAT * THE ELEMENTS OF THE RIGHT-HAND SIDE VECTOR Y ARE VERY * LARGE. * * ===================================================================== * * .. PARAMETERS .. DOUBLE PRECISION ONE, ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. LOCAL SCALARS .. INTEGER K DOUBLE PRECISION ABSAK, AK, BIGNUM, EPS, PERT, SFMIN, TEMP * .. * .. INTRINSIC FUNCTIONS .. INTRINSIC ABS, MAX, SIGN * .. * .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION DLAMCH EXTERNAL DLAMCH * .. * .. EXTERNAL SUBROUTINES .. EXTERNAL XERBLA * .. * .. EXECUTABLE STATEMENTS .. * INFO = 0 IF( ( ABS( JOB ).GT.2 ) .OR. ( JOB.EQ.0 ) ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DLAGTS', -INFO ) RETURN END IF * IF( N.EQ.0 ) $ RETURN * EPS = DLAMCH( 'EPSILON' ) SFMIN = DLAMCH( 'SAFE MINIMUM' ) BIGNUM = ONE / SFMIN * IF( JOB.LT.0 ) THEN IF( TOL.LE.ZERO ) THEN TOL = ABS( A( 1 ) ) IF( N.GT.1 ) $ TOL = MAX( TOL, ABS( A( 2 ) ), ABS( B( 1 ) ) ) DO 10 K = 3, N TOL = MAX( TOL, ABS( A( K ) ), ABS( B( K-1 ) ), $ ABS( D( K-2 ) ) ) 10 CONTINUE TOL = TOL*EPS IF( TOL.EQ.ZERO ) $ TOL = EPS END IF END IF * IF( ABS( JOB ).EQ.1 ) THEN DO 20 K = 2, N IF( IN( K-1 ).EQ.0 ) THEN Y( K ) = Y( K ) - C( K-1 )*Y( K-1 ) ELSE TEMP = Y( K-1 ) Y( K-1 ) = Y( K ) Y( K ) = TEMP - C( K-1 )*Y( K ) END IF 20 CONTINUE IF( JOB.EQ.1 ) THEN DO 30 K = N, 1, -1 IF( K.LE.N-2 ) THEN TEMP = Y( K ) - B( K )*Y( K+1 ) - D( K )*Y( K+2 ) ELSE IF( K.EQ.N-1 ) THEN TEMP = Y( K ) - B( K )*Y( K+1 ) ELSE TEMP = Y( K ) END IF AK = A( K ) ABSAK = ABS( AK ) IF( ABSAK.LT.ONE ) THEN IF( ABSAK.LT.SFMIN ) THEN IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK ) $ THEN INFO = K RETURN ELSE TEMP = TEMP*BIGNUM AK = AK*BIGNUM END IF ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN INFO = K RETURN END IF END IF Y( K ) = TEMP / AK 30 CONTINUE ELSE DO 50 K = N, 1, -1 IF( K.LE.N-2 ) THEN TEMP = Y( K ) - B( K )*Y( K+1 ) - D( K )*Y( K+2 ) ELSE IF( K.EQ.N-1 ) THEN TEMP = Y( K ) - B( K )*Y( K+1 ) ELSE TEMP = Y( K ) END IF AK = A( K ) PERT = SIGN( TOL, AK ) 40 CONTINUE ABSAK = ABS( AK ) IF( ABSAK.LT.ONE ) THEN IF( ABSAK.LT.SFMIN ) THEN IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK ) $ THEN AK = AK + PERT PERT = 2*PERT GO TO 40 ELSE TEMP = TEMP*BIGNUM AK = AK*BIGNUM END IF ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN AK = AK + PERT PERT = 2*PERT GO TO 40 END IF END IF Y( K ) = TEMP / AK 50 CONTINUE END IF ELSE * * COME TO HERE IF JOB = 2 OR -2 * IF( JOB.EQ.2 ) THEN DO 60 K = 1, N IF( K.GE.3 ) THEN TEMP = Y( K ) - B( K-1 )*Y( K-1 ) - D( K-2 )*Y( K-2 ) ELSE IF( K.EQ.2 ) THEN TEMP = Y( K ) - B( K-1 )*Y( K-1 ) ELSE TEMP = Y( K ) END IF AK = A( K ) ABSAK = ABS( AK ) IF( ABSAK.LT.ONE ) THEN IF( ABSAK.LT.SFMIN ) THEN IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK ) $ THEN INFO = K RETURN ELSE TEMP = TEMP*BIGNUM AK = AK*BIGNUM END IF ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN INFO = K RETURN END IF END IF Y( K ) = TEMP / AK 60 CONTINUE ELSE DO 80 K = 1, N IF( K.GE.3 ) THEN TEMP = Y( K ) - B( K-1 )*Y( K-1 ) - D( K-2 )*Y( K-2 ) ELSE IF( K.EQ.2 ) THEN TEMP = Y( K ) - B( K-1 )*Y( K-1 ) ELSE TEMP = Y( K ) END IF AK = A( K ) PERT = SIGN( TOL, AK ) 70 CONTINUE ABSAK = ABS( AK ) IF( ABSAK.LT.ONE ) THEN IF( ABSAK.LT.SFMIN ) THEN IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK ) $ THEN AK = AK + PERT PERT = 2*PERT GO TO 70 ELSE TEMP = TEMP*BIGNUM AK = AK*BIGNUM END IF ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN AK = AK + PERT PERT = 2*PERT GO TO 70 END IF END IF Y( K ) = TEMP / AK 80 CONTINUE END IF * DO 90 K = N, 2, -1 IF( IN( K-1 ).EQ.0 ) THEN Y( K-1 ) = Y( K-1 ) - C( K-1 )*Y( K ) ELSE TEMP = Y( K-1 ) Y( K-1 ) = Y( K ) Y( K ) = TEMP - C( K-1 )*Y( K ) END IF 90 CONTINUE END IF * * END OF DLAGTS * END CUT HERE............ CAT > DLAGTF.F <<'CUT HERE............' SUBROUTINE DLAGTF( N, A, LAMBDA, B, C, TOL, D, IN, INFO ) * * -- LAPACK ROUTINE (VERSION 1.1) -- * UNIV. OF TENNESSEE, UNIV. OF CALIFORNIA BERKELEY, NAG LTD., * COURANT INSTITUTE, ARGONNE NATIONAL LAB, AND RICE UNIVERSITY * OCTOBER 31, 1992 * * .. SCALAR ARGUMENTS .. INTEGER INFO, N DOUBLE PRECISION LAMBDA, TOL * .. * .. ARRAY ARGUMENTS .. INTEGER IN( * ) DOUBLE PRECISION A( * ), B( * ), C( * ), D( * ) * .. * * PURPOSE * ======= * * DLAGTF FACTORIZES THE MATRIX (T - LAMBDA*I), WHERE T IS AN N BY N * TRIDIAGONAL MATRIX AND LAMBDA IS A SCALAR, AS * * T - LAMBDA*I = PLU, * * WHERE P IS A PERMUTATION MATRIX, L IS A UNIT LOWER TRIDIAGONAL MATRIX * WITH AT MOST ONE NON-ZERO SUB-DIAGONAL ELEMENTS PER COLUMN AND U IS * AN UPPER TRIANGULAR MATRIX WITH AT MOST TWO NON-ZERO SUPER-DIAGONAL * ELEMENTS PER COLUMN. * * THE FACTORIZATION IS OBTAINED BY GAUSSIAN ELIMINATION WITH PARTIAL * PIVOTING AND IMPLICIT ROW SCALING. * * THE PARAMETER LAMBDA IS INCLUDED IN THE ROUTINE SO THAT DLAGTF MAY * BE USED, IN CONJUNCTION WITH DLAGTS, TO OBTAIN EIGENVECTORS OF T BY * INVERSE ITERATION. * * ARGUMENTS * ========= * * N (INPUT) INTEGER * THE ORDER OF THE MATRIX T. * * A (INPUT/OUTPUT) DOUBLE PRECISION ARRAY, DIMENSION (N) * ON ENTRY, A MUST CONTAIN THE DIAGONAL ELEMENTS OF T. * * ON EXIT, A IS OVERWRITTEN BY THE N DIAGONAL ELEMENTS OF THE * UPPER TRIANGULAR MATRIX U OF THE FACTORIZATION OF T. * * LAMBDA (INPUT) DOUBLE PRECISION * ON ENTRY, THE SCALAR LAMBDA. * * B (INPUT/OUTPUT) DOUBLE PRECISION ARRAY, DIMENSION (N-1) * ON ENTRY, B MUST CONTAIN THE (N-1) SUPER-DIAGONAL ELEMENTS OF * T. * * ON EXIT, B IS OVERWRITTEN BY THE (N-1) SUPER-DIAGONAL * ELEMENTS OF THE MATRIX U OF THE FACTORIZATION OF T. * * C (INPUT/OUTPUT) DOUBLE PRECISION ARRAY, DIMENSION (N-1) * ON ENTRY, C MUST CONTAIN THE (N-1) SUB-DIAGONAL ELEMENTS OF * T. * * ON EXIT, C IS OVERWRITTEN BY THE (N-1) SUB-DIAGONAL ELEMENTS * OF THE MATRIX L OF THE FACTORIZATION OF T. * * TOL (INPUT) DOUBLE PRECISION * ON ENTRY, A RELATIVE TOLERANCE USED TO INDICATE WHETHER OR * NOT THE MATRIX (T - LAMBDA*I) IS NEARLY SINGULAR. TOL SHOULD * NORMALLY BE CHOSE AS APPROXIMATELY THE LARGEST RELATIVE ERROR * IN THE ELEMENTS OF T. FOR EXAMPLE, IF THE ELEMENTS OF T ARE * CORRECT TO ABOUT 4 SIGNIFICANT FIGURES, THEN TOL SHOULD BE * SET TO ABOUT 5*10**(-4). IF TOL IS SUPPLIED AS LESS THAN EPS, * WHERE EPS IS THE RELATIVE MACHINE PRECISION, THEN THE VALUE * EPS IS USED IN PLACE OF TOL. * * D (OUTPUT) DOUBLE PRECISION ARRAY, DIMENSION (N-2) * ON EXIT, D IS OVERWRITTEN BY THE (N-2) SECOND SUPER-DIAGONAL * ELEMENTS OF THE MATRIX U OF THE FACTORIZATION OF T. * * IN (OUTPUT) INTEGER ARRAY, DIMENSION (N) * ON EXIT, IN CONTAINS DETAILS OF THE PERMUTATION MATRIX P. IF * AN INTERCHANGE OCCURRED AT THE KTH STEP OF THE ELIMINATION, * THEN IN(K) = 1, OTHERWISE IN(K) = 0. THE ELEMENT IN(N) * RETURNS THE SMALLEST POSITIVE INTEGER J SUCH THAT * * ABS( U(J,J) ).LE. NORM( (T - LAMBDA*I)(J) )*TOL, * * WHERE NORM( A(J) ) DENOTES THE SUM OF THE ABSOLUTE VALUES OF * THE JTH ROW OF THE MATRIX A. IF NO SUCH J EXISTS THEN IN(N) * IS RETURNED AS ZERO. IF IN(N) IS RETURNED AS POSITIVE, THEN A * DIAGONAL ELEMENT OF U IS SMALL, INDICATING THAT * (T - LAMBDA*I) IS SINGULAR OR NEARLY SINGULAR, * * INFO (OUTPUT) * = 0 : SUCCESSFUL EXIT * .LT. 0: IF INFO = -K, THE KTH ARGUMENT HAD AN ILLEGAL VALUE * * ===================================================================== * * .. PARAMETERS .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) * .. * .. LOCAL SCALARS .. INTEGER K DOUBLE PRECISION EPS, MULT, PIV1, PIV2, SCALE1, SCALE2, TEMP, TL * .. * .. INTRINSIC FUNCTIONS .. INTRINSIC ABS, MAX * .. * .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION DLAMCH EXTERNAL DLAMCH * .. * .. EXTERNAL SUBROUTINES .. EXTERNAL XERBLA * .. * .. EXECUTABLE STATEMENTS .. * INFO = 0 IF( N.LT.0 ) THEN INFO = -1 CALL XERBLA( 'DLAGTF', -INFO ) RETURN END IF * IF( N.EQ.0 ) $ RETURN * A( 1 ) = A( 1 ) - LAMBDA IN( N ) = 0 IF( N.EQ.1 ) THEN IF( A( 1 ).EQ.ZERO ) $ IN( 1 ) = 1 RETURN END IF * EPS = DLAMCH( 'EPSILON' ) * TL = MAX( TOL, EPS ) SCALE1 = ABS( A( 1 ) ) + ABS( B( 1 ) ) DO 10 K = 1, N - 1 A( K+1 ) = A( K+1 ) - LAMBDA SCALE2 = ABS( C( K ) ) + ABS( A( K+1 ) ) IF( K.LT.( N-1 ) ) $ SCALE2 = SCALE2 + ABS( B( K+1 ) ) IF( A( K ).EQ.ZERO ) THEN PIV1 = ZERO ELSE PIV1 = ABS( A( K ) ) / SCALE1 END IF IF( C( K ).EQ.ZERO ) THEN IN( K ) = 0 PIV2 = ZERO SCALE1 = SCALE2 IF( K.LT.( N-1 ) ) $ D( K ) = ZERO ELSE PIV2 = ABS( C( K ) ) / SCALE2 IF( PIV2.LE.PIV1 ) THEN IN( K ) = 0 SCALE1 = SCALE2 C( K ) = C( K ) / A( K ) A( K+1 ) = A( K+1 ) - C( K )*B( K ) IF( K.LT.( N-1 ) ) $ D( K ) = ZERO ELSE IN( K ) = 1 MULT = A( K ) / C( K ) A( K ) = C( K ) TEMP = A( K+1 ) A( K+1 ) = B( K ) - MULT*TEMP IF( K.LT.( N-1 ) ) THEN D( K ) = B( K+1 ) B( K+1 ) = -MULT*D( K ) END IF B( K ) = TEMP C( K ) = MULT END IF END IF IF( ( MAX( PIV1, PIV2 ).LE.TL ) .AND. ( IN( N ).EQ.0 ) ) $ IN( N ) = K 10 CONTINUE IF( ( ABS( A( N ) ).LE.SCALE1*TL ) .AND. ( IN( N ).EQ.0 ) ) $ IN( N ) = N * RETURN * * END OF DLAGTF * END CUT HERE............ CAT > DLARNV.F <<'CUT HERE............' SUBROUTINE DLARNV( IDIST, ISEED, N, X ) * * -- LAPACK AUXILIARY ROUTINE (VERSION 1.1) -- * UNIV. OF TENNESSEE, UNIV. OF CALIFORNIA BERKELEY, NAG LTD., * COURANT INSTITUTE, ARGONNE NATIONAL LAB, AND RICE UNIVERSITY * OCTOBER 31, 1992 * * .. SCALAR ARGUMENTS .. INTEGER IDIST, N * .. * .. ARRAY ARGUMENTS .. INTEGER ISEED( 4 ) DOUBLE PRECISION X( * ) * .. * * PURPOSE * ======= * * DLARNV RETURNS A VECTOR OF N RANDOM REAL NUMBERS FROM A UNIFORM OR * NORMAL DISTRIBUTION. * * ARGUMENTS * ========= * * IDIST (INPUT) INTEGER * SPECIFIES THE DISTRIBUTION OF THE RANDOM NUMBERS: * = 1: UNIFORM (0,1) * = 2: UNIFORM (-1,1) * = 3: NORMAL (0,1) * * ISEED (INPUT/OUTPUT) INTEGER ARRAY, DIMENSION (4) * ON ENTRY, THE SEED OF THE RANDOM NUMBER GENERATOR; THE ARRAY * ELEMENTS MUST BE BETWEEN 0 AND 4095, AND ISEED(4) MUST BE * ODD. * ON EXIT, THE SEED IS UPDATED. * * N (INPUT) INTEGER * THE NUMBER OF RANDOM NUMBERS TO BE GENERATED. * * X (OUTPUT) DOUBLE PRECISION ARRAY, DIMENSION (N) * THE GENERATED RANDOM NUMBERS. * * FURTHER DETAILS * =============== * * THIS ROUTINE CALLS THE AUXILIARY ROUTINE DLARUV TO GENERATE RANDOM * REAL NUMBERS FROM A UNIFORM (0,1) DISTRIBUTION, IN BATCHES OF UP TO * 128 USING VECTORISABLE CODE. THE BOX-MULLER METHOD IS USED TO * TRANSFORM NUMBERS FROM A UNIFORM TO A NORMAL DISTRIBUTION. * * ===================================================================== * * .. PARAMETERS .. DOUBLE PRECISION ONE, TWO PARAMETER ( ONE = 1.0D+0, TWO = 2.0D+0 ) INTEGER LV PARAMETER ( LV = 128 ) DOUBLE PRECISION TWOPI PARAMETER ( TWOPI = 6.28318530717958623199592D+0 ) * .. * .. LOCAL SCALARS .. INTEGER I, IL, IL2, IV * .. * .. LOCAL ARRAYS .. DOUBLE PRECISION U( LV ) * .. * .. INTRINSIC FUNCTIONS .. INTRINSIC COS, LOG, MIN, SQRT * .. * .. EXTERNAL SUBROUTINES .. EXTERNAL DLARUV * .. * .. EXECUTABLE STATEMENTS .. * DO 40 IV = 1, N, LV / 2 IL = MIN( LV / 2, N-IV+1 ) IF( IDIST.EQ.3 ) THEN IL2 = 2*IL ELSE IL2 = IL END IF * * CALL DLARUV TO GENERATE IL2 NUMBERS FROM A UNIFORM (0,1) * DISTRIBUTION (IL2 <= LV) * CALL DLARUV( ISEED, IL2, U ) * IF( IDIST.EQ.1 ) THEN * * COPY GENERATED NUMBERS * DO 10 I = 1, IL X( IV+I-1 ) = U( I ) 10 CONTINUE ELSE IF( IDIST.EQ.2 ) THEN * * CONVERT GENERATED NUMBERS TO UNIFORM (-1,1) DISTRIBUTION * DO 20 I = 1, IL X( IV+I-1 ) = TWO*U( I ) - ONE 20 CONTINUE ELSE IF( IDIST.EQ.3 ) THEN * * CONVERT GENERATED NUMBERS TO NORMAL (0,1) DISTRIBUTION * DO 30 I = 1, IL X( IV+I-1 ) = SQRT( -TWO*LOG( U( 2*I-1 ) ) )* $ COS( TWOPI*U( 2*I ) ) 30 CONTINUE END IF 40 CONTINUE RETURN * * END OF DLARNV * END CUT HERE............ CAT > DLARUV.F <<'CUT HERE............' SUBROUTINE DLARUV( ISEED, N, X ) * * -- LAPACK AUXILIARY ROUTINE (VERSION 1.1) -- * UNIV. OF TENNESSEE, UNIV. OF CALIFORNIA BERKELEY, NAG LTD., * COURANT INSTITUTE, ARGONNE NATIONAL LAB, AND RICE UNIVERSITY * OCTOBER 31, 1992 * * .. SCALAR ARGUMENTS .. INTEGER N * .. * .. ARRAY ARGUMENTS .. INTEGER ISEED( 4 ) DOUBLE PRECISION X( N ) * .. * * PURPOSE * ======= * * DLARUV RETURNS A VECTOR OF N RANDOM REAL NUMBERS FROM A UNIFORM (0,1) * DISTRIBUTION (N <= 128). * * THIS IS AN AUXILIARY ROUTINE CALLED BY DLARNV AND ZLARNV. * * ARGUMENTS * ========= * * ISEED (INPUT/OUTPUT) INTEGER ARRAY, DIMENSION (4) * ON ENTRY, THE SEED OF THE RANDOM NUMBER GENERATOR; THE ARRAY * ELEMENTS MUST BE BETWEEN 0 AND 4095, AND ISEED(4) MUST BE * ODD. * ON EXIT, THE SEED IS UPDATED. * * N (INPUT) INTEGER * THE NUMBER OF RANDOM NUMBERS TO BE GENERATED. N <= 128. * * X (OUTPUT) DOUBLE PRECISION ARRAY, DIMENSION (N) * THE GENERATED RANDOM NUMBERS. * * FURTHER DETAILS * =============== * * THIS ROUTINE USES A MULTIPLICATIVE CONGRUENTIAL METHOD WITH MODULUS * 2**48 AND MULTIPLIER 33952834046453 (SEE G.S.FISHMAN, * 'MULTIPLICATIVE CONGRUENTIAL RANDOM NUMBER GENERATORS WITH MODULUS * 2**B: AN EXHAUSTIVE ANALYSIS FOR B = 32 AND A PARTIAL ANALYSIS FOR * B = 48', MATH. COMP. 189, PP 331-344, 1990). * * 48-BIT INTEGERS ARE STORED IN 4 INTEGER ARRAY ELEMENTS WITH 12 BITS * PER ELEMENT. HENCE THE ROUTINE IS PORTABLE ACROSS MACHINES WITH * INTEGERS OF 32 BITS OR MORE. * * ===================================================================== * * .. PARAMETERS .. DOUBLE PRECISION ONE PARAMETER ( ONE = 1.0D0 ) INTEGER LV, IPW2 DOUBLE PRECISION R PARAMETER ( LV = 128, IPW2 = 4096, R = ONE / IPW2 ) * .. * .. LOCAL SCALARS .. INTEGER I, I1, I2, I3, I4, IT1, IT2, IT3, IT4, J * .. * .. LOCAL ARRAYS .. INTEGER MM( LV, 4 ) * .. * .. INTRINSIC FUNCTIONS .. INTRINSIC DBLE, MIN, MOD * .. * .. DATA STATEMENTS .. DATA ( MM( 1, J ), J = 1, 4 ) / 494, 322, 2508, $ 2549 / DATA ( MM( 2, J ), J = 1, 4 ) / 2637, 789, 3754, $ 1145 / DATA ( MM( 3, J ), J = 1, 4 ) / 255, 1440, 1766, $ 2253 / DATA ( MM( 4, J ), J = 1, 4 ) / 2008, 752, 3572, $ 305 / DATA ( MM( 5, J ), J = 1, 4 ) / 1253, 2859, 2893, $ 3301 / DATA ( MM( 6, J ), J = 1, 4 ) / 3344, 123, 307, $ 1065 / DATA ( MM( 7, J ), J = 1, 4 ) / 4084, 1848, 1297, $ 3133 / DATA ( MM( 8, J ), J = 1, 4 ) / 1739, 643, 3966, $ 2913 / DATA ( MM( 9, J ), J = 1, 4 ) / 3143, 2405, 758, $ 3285 / DATA ( MM( 10, J ), J = 1, 4 ) / 3468, 2638, 2598, $ 1241 / DATA ( MM( 11, J ), J = 1, 4 ) / 688, 2344, 3406, $ 1197 / DATA ( MM( 12, J ), J = 1, 4 ) / 1657, 46, 2922, $ 3729 / DATA ( MM( 13, J ), J = 1, 4 ) / 1238, 3814, 1038, $ 2501 / DATA ( MM( 14, J ), J = 1, 4 ) / 3166, 913, 2934, $ 1673 / DATA ( MM( 15, J ), J = 1, 4 ) / 1292, 3649, 2091, $ 541 / DATA ( MM( 16, J ), J = 1, 4 ) / 3422, 339, 2451, $ 2753 / DATA ( MM( 17, J ), J = 1, 4 ) / 1270, 3808, 1580, $ 949 / DATA ( MM( 18, J ), J = 1, 4 ) / 2016, 822, 1958, $ 2361 / DATA ( MM( 19, J ), J = 1, 4 ) / 154, 2832, 2055, $ 1165 / DATA ( MM( 20, J ), J = 1, 4 ) / 2862, 3078, 1507, $ 4081 / DATA ( MM( 21, J ), J = 1, 4 ) / 697, 3633, 1078, $ 2725 / DATA ( MM( 22, J ), J = 1, 4 ) / 1706, 2970, 3273, $ 3305 / DATA ( MM( 23, J ), J = 1, 4 ) / 491, 637, 17, $ 3069 / DATA ( MM( 24, J ), J = 1, 4 ) / 931, 2249, 854, $ 3617 / DATA ( MM( 25, J ), J = 1, 4 ) / 1444, 2081, 2916, $ 3733 / DATA ( MM( 26, J ), J = 1, 4 ) / 444, 4019, 3971, $ 409 / DATA ( MM( 27, J ), J = 1, 4 ) / 3577, 1478, 2889, $ 2157 / DATA ( MM( 28, J ), J = 1, 4 ) / 3944, 242, 3831, $ 1361 / DATA ( MM( 29, J ), J = 1, 4 ) / 2184, 481, 2621, $ 3973 / DATA ( MM( 30, J ), J = 1, 4 ) / 1661, 2075, 1541, $ 1865 / DATA ( MM( 31, J ), J = 1, 4 ) / 3482, 4058, 893, $ 2525 / DATA ( MM( 32, J ), J = 1, 4 ) / 657, 622, 736, $ 1409 / DATA ( MM( 33, J ), J = 1, 4 ) / 3023, 3376, 3992, $ 3445 / DATA ( MM( 34, J ), J = 1, 4 ) / 3618, 812, 787, $ 3577 / DATA ( MM( 35, J ), J = 1, 4 ) / 1267, 234, 2125, $ 77 / DATA ( MM( 36, J ), J = 1, 4 ) / 1828, 641, 2364, $ 3761 / DATA ( MM( 37, J ), J = 1, 4 ) / 164, 4005, 2460, $ 2149 / DATA ( MM( 38, J ), J = 1, 4 ) / 3798, 1122, 257, $ 1449 / DATA ( MM( 39, J ), J = 1, 4 ) / 3087, 3135, 1574, $ 3005 / DATA ( MM( 40, J ), J = 1, 4 ) / 2400, 2640, 3912, $ 225 / DATA ( MM( 41, J ), J = 1, 4 ) / 2870, 2302, 1216, $ 85 / DATA ( MM( 42, J ), J = 1, 4 ) / 3876, 40, 3248, $ 3673 / DATA ( MM( 43, J ), J = 1, 4 ) / 1905, 1832, 3401, $ 3117 / DATA ( MM( 44, J ), J = 1, 4 ) / 1593, 2247, 2124, $ 3089 / DATA ( MM( 45, J ), J = 1, 4 ) / 1797, 2034, 2762, $ 1349 / DATA ( MM( 46, J ), J = 1, 4 ) / 1234, 2637, 149, $ 2057 / DATA ( MM( 47, J ), J = 1, 4 ) / 3460, 1287, 2245, $ 413 / DATA ( MM( 48, J ), J = 1, 4 ) / 328, 1691, 166, $ 65 / DATA ( MM( 49, J ), J = 1, 4 ) / 2861, 496, 466, $ 1845 / DATA ( MM( 50, J ), J = 1, 4 ) / 1950, 1597, 4018, $ 697 / DATA ( MM( 51, J ), J = 1, 4 ) / 617, 2394, 1399, $ 3085 / DATA ( MM( 52, J ), J = 1, 4 ) / 2070, 2584, 190, $ 3441 / DATA ( MM( 53, J ), J = 1, 4 ) / 3331, 1843, 2879, $ 1573 / DATA ( MM( 54, J ), J = 1, 4 ) / 769, 336, 153, $ 3689 / DATA ( MM( 55, J ), J = 1, 4 ) / 1558, 1472, 2320, $ 2941 / DATA ( MM( 56, J ), J = 1, 4 ) / 2412, 2407, 18, $ 929 / DATA ( MM( 57, J ), J = 1, 4 ) / 2800, 433, 712, $ 533 / DATA ( MM( 58, J ), J = 1, 4 ) / 189, 2096, 2159, $ 2841 / DATA ( MM( 59, J ), J = 1, 4 ) / 287, 1761, 2318, $ 4077 / DATA ( MM( 60, J ), J = 1, 4 ) / 2045, 2810, 2091, $ 721 / DATA ( MM( 61, J ), J = 1, 4 ) / 1227, 566, 3443, $ 2821 / DATA ( MM( 62, J ), J = 1, 4 ) / 2838, 442, 1510, $ 2249 / DATA ( MM( 63, J ), J = 1, 4 ) / 209, 41, 449, $ 2397 / DATA ( MM( 64, J ), J = 1, 4 ) / 2770, 1238, 1956, $ 2817 / DATA ( MM( 65, J ), J = 1, 4 ) / 3654, 1086, 2201, $ 245 / DATA ( MM( 66, J ), J = 1, 4 ) / 3993, 603, 3137, $ 1913 / DATA ( MM( 67, J ), J = 1, 4 ) / 192, 840, 3399, $ 1997 / DATA ( MM( 68, J ), J = 1, 4 ) / 2253, 3168, 1321, $ 3121 / DATA ( MM( 69, J ), J = 1, 4 ) / 3491, 1499, 2271, $ 997 / DATA ( MM( 70, J ), J = 1, 4 ) / 2889, 1084, 3667, $ 1833 / DATA ( MM( 71, J ), J = 1, 4 ) / 2857, 3438, 2703, $ 2877 / DATA ( MM( 72, J ), J = 1, 4 ) / 2094, 2408, 629, $ 1633 / DATA ( MM( 73, J ), J = 1, 4 ) / 1818, 1589, 2365, $ 981 / DATA ( MM( 74, J ), J = 1, 4 ) / 688, 2391, 2431, $ 2009 / DATA ( MM( 75, J ), J = 1, 4 ) / 1407, 288, 1113, $ 941 / DATA ( MM( 76, J ), J = 1, 4 ) / 634, 26, 3922, $ 2449 / DATA ( MM( 77, J ), J = 1, 4 ) / 3231, 512, 2554, $ 197 / DATA ( MM( 78, J ), J = 1, 4 ) / 815, 1456, 184, $ 2441 / DATA ( MM( 79, J ), J = 1, 4 ) / 3524, 171, 2099, $ 285 / DATA ( MM( 80, J ), J = 1, 4 ) / 1914, 1677, 3228, $ 1473 / DATA ( MM( 81, J ), J = 1, 4 ) / 516, 2657, 4012, $ 2741 / DATA ( MM( 82, J ), J = 1, 4 ) / 164, 2270, 1921, $ 3129 / DATA ( MM( 83, J ), J = 1, 4 ) / 303, 2587, 3452, $ 909 / DATA ( MM( 84, J ), J = 1, 4 ) / 2144, 2961, 3901, $ 2801 / DATA ( MM( 85, J ), J = 1, 4 ) / 3480, 1970, 572, $ 421 / DATA ( MM( 86, J ), J = 1, 4 ) / 119, 1817, 3309, $ 4073 / DATA ( MM( 87, J ), J = 1, 4 ) / 3357, 676, 3171, $ 2813 / DATA ( MM( 88, J ), J = 1, 4 ) / 837, 1410, 817, $ 2337 / DATA ( MM( 89, J ), J = 1, 4 ) / 2826, 3723, 3039, $ 1429 / DATA ( MM( 90, J ), J = 1, 4 ) / 2332, 2803, 1696, $ 1177 / DATA ( MM( 91, J ), J = 1, 4 ) / 2089, 3185, 1256, $ 1901 / DATA ( MM( 92, J ), J = 1, 4 ) / 3780, 184, 3715, $ 81 / DATA ( MM( 93, J ), J = 1, 4 ) / 1700, 663, 2077, $ 1669 / DATA ( MM( 94, J ), J = 1, 4 ) / 3712, 499, 3019, $ 2633 / DATA ( MM( 95, J ), J = 1, 4 ) / 150, 3784, 1497, $ 2269 / DATA ( MM( 96, J ), J = 1, 4 ) / 2000, 1631, 1101, $ 129 / DATA ( MM( 97, J ), J = 1, 4 ) / 3375, 1925, 717, $ 1141 / DATA ( MM( 98, J ), J = 1, 4 ) / 1621, 3912, 51, $ 249 / DATA ( MM( 99, J ), J = 1, 4 ) / 3090, 1398, 981, $ 3917 / DATA ( MM( 100, J ), J = 1, 4 ) / 3765, 1349, 1978, $ 2481 / DATA ( MM( 101, J ), J = 1, 4 ) / 1149, 1441, 1813, $ 3941 / DATA ( MM( 102, J ), J = 1, 4 ) / 3146, 2224, 3881, $ 2217 / DATA ( MM( 103, J ), J = 1, 4 ) / 33, 2411, 76, $ 2749 / DATA ( MM( 104, J ), J = 1, 4 ) / 3082, 1907, 3846, $ 3041 / DATA ( MM( 105, J ), J = 1, 4 ) / 2741, 3192, 3694, $ 1877 / DATA ( MM( 106, J ), J = 1, 4 ) / 359, 2786, 1682, $ 345 / DATA ( MM( 107, J ), J = 1, 4 ) / 3316, 382, 124, $ 2861 / DATA ( MM( 108, J ), J = 1, 4 ) / 1749, 37, 1660, $ 1809 / DATA ( MM( 109, J ), J = 1, 4 ) / 185, 759, 3997, $ 3141 / DATA ( MM( 110, J ), J = 1, 4 ) / 2784, 2948, 479, $ 2825 / DATA ( MM( 111, J ), J = 1, 4 ) / 2202, 1862, 1141, $ 157 / DATA ( MM( 112, J ), J = 1, 4 ) / 2199, 3802, 886, $ 2881 / DATA ( MM( 113, J ), J = 1, 4 ) / 1364, 2423, 3514, $ 3637 / DATA ( MM( 114, J ), J = 1, 4 ) / 1244, 2051, 1301, $ 1465 / DATA ( MM( 115, J ), J = 1, 4 ) / 2020, 2295, 3604, $ 2829 / DATA ( MM( 116, J ), J = 1, 4 ) / 3160, 1332, 1888, $ 2161 / DATA ( MM( 117, J ), J = 1, 4 ) / 2785, 1832, 1836, $ 3365 / DATA ( MM( 118, J ), J = 1, 4 ) / 2772, 2405, 1990, $ 361 / DATA ( MM( 119, J ), J = 1, 4 ) / 1217, 3638, 2058, $ 2685 / DATA ( MM( 120, J ), J = 1, 4 ) / 1822, 3661, 692, $ 3745 / DATA ( MM( 121, J ), J = 1, 4 ) / 1245, 327, 1194, $ 2325 / DATA ( MM( 122, J ), J = 1, 4 ) / 2252, 3660, 20, $ 3609 / DATA ( MM( 123, J ), J = 1, 4 ) / 3904, 716, 3285, $ 3821 / DATA ( MM( 124, J ), J = 1, 4 ) / 2774, 1842, 2046, $ 3537 / DATA ( MM( 125, J ), J = 1, 4 ) / 997, 3987, 2107, $ 517 / DATA ( MM( 126, J ), J = 1, 4 ) / 2573, 1368, 3508, $ 3017 / DATA ( MM( 127, J ), J = 1, 4 ) / 1148, 1848, 3525, $ 2141 / DATA ( MM( 128, J ), J = 1, 4 ) / 545, 2366, 3801, $ 1537 / * .. * .. EXECUTABLE STATEMENTS .. * I1 = ISEED( 1 ) I2 = ISEED( 2 ) I3 = ISEED( 3 ) I4 = ISEED( 4 ) * DO 10 I = 1, MIN( N, LV ) * * MULTIPLY THE SEED BY I-TH POWER OF THE MULTIPLIER MODULO 2**48 * IT4 = I4*MM( I, 4 ) IT3 = IT4 / IPW2 IT4 = IT4 - IPW2*IT3 IT3 = IT3 + I3*MM( I, 4 ) + I4*MM( I, 3 ) IT2 = IT3 / IPW2 IT3 = IT3 - IPW2*IT2 IT2 = IT2 + I2*MM( I, 4 ) + I3*MM( I, 3 ) + I4*MM( I, 2 ) IT1 = IT2 / IPW2 IT2 = IT2 - IPW2*IT1 IT1 = IT1 + I1*MM( I, 4 ) + I2*MM( I, 3 ) + I3*MM( I, 2 ) + $ I4*MM( I, 1 ) IT1 = MOD( IT1, IPW2 ) * * CONVERT 48-BIT INTEGER TO A REAL NUMBER IN THE INTERVAL (0,1) * X( I ) = R*( DBLE( IT1 )+R*( DBLE( IT2 )+R*( DBLE( IT3 )+R* $ DBLE( IT4 ) ) ) ) 10 CONTINUE * * RETURN FINAL VALUE OF SEED * ISEED( 1 ) = IT1 ISEED( 2 ) = IT2 ISEED( 3 ) = IT3 ISEED( 4 ) = IT4 RETURN * * END OF DLARUV * END CUT HERE............ CAT > DSTEBZ.F <<'CUT HERE............' SUBROUTINE DSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, $ M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, $ INFO ) * * -- LAPACK ROUTINE (VERSION 1.1) -- * UNIV. OF TENNESSEE, UNIV. OF CALIFORNIA BERKELEY, NAG LTD., * COURANT INSTITUTE, ARGONNE NATIONAL LAB, AND RICE UNIVERSITY * MARCH 31, 1993 * * .. SCALAR ARGUMENTS .. CHARACTER ORDER, RANGE INTEGER IL, INFO, IU, M, N, NSPLIT DOUBLE PRECISION ABSTOL, VL, VU * .. * .. ARRAY ARGUMENTS .. INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * ) DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ) * .. * * PURPOSE * ======= * * DSTEBZ COMPUTES THE EIGENVALUES OF A SYMMETRIC TRIDIAGONAL * MATRIX T. THE USER MAY ASK FOR ALL EIGENVALUES, ALL EIGENVALUES * IN THE HALF-OPEN INTERVAL (VL, VU], OR THE IL-TH THROUGH IU-TH * EIGENVALUES. * * SEE W. KAHAN "ACCURATE EIGENVALUES OF A SYMMETRIC TRIDIAGONAL * MATRIX", REPORT CS41, COMPUTER SCIENCE DEPT., STANFORD * UNIVERSITY, JULY 21, 1966. * * ARGUMENTS * ========= * * RANGE (INPUT) CHARACTER * = 'A': ("ALL") ALL EIGENVALUES WILL BE FOUND. * = 'V': ("VALUE") ALL EIGENVALUES IN THE HALF-OPEN INTERVAL * (VL, VU] WILL BE FOUND. * = 'I': ("INDEX") THE IL-TH THROUGH IU-TH EIGENVALUES (OF THE * ENTIRE MATRIX) WILL BE FOUND. * * ORDER (INPUT) CHARACTER * = 'B': ("BY BLOCK") THE EIGENVALUES WILL BE GROUPED BY * SPLIT-OFF BLOCK (SEE IBLOCK, ISPLIT) AND * ORDERED FROM SMALLEST TO LARGEST WITHIN * THE BLOCK. * = 'E': ("ENTIRE MATRIX") * THE EIGENVALUES FOR THE ENTIRE MATRIX * WILL BE ORDERED FROM SMALLEST TO * LARGEST. * * N (INPUT) INTEGER * THE DIMENSION OF THE TRIDIAGONAL MATRIX T. N >= 0. * * VL (INPUT) DOUBLE PRECISION * IF RANGE='V', THE LOWER BOUND OF THE INTERVAL TO BE SEARCHED * FOR EIGENVALUES. EIGENVALUES LESS THAN OR EQUAL TO VL WILL * NOT BE RETURNED. NOT REFERENCED IF RANGE='A' OR 'I'. * * VU (INPUT) DOUBLE PRECISION * IF RANGE='V', THE UPPER BOUND OF THE INTERVAL TO BE SEARCHED * FOR EIGENVALUES. EIGENVALUES GREATER THAN VU WILL NOT BE * RETURNED. VU MUST BE GREATER THAN VL. NOT REFERENCED IF * RANGE='A' OR 'I'. * * IL (INPUT) INTEGER * IF RANGE='I', THE INDEX (FROM SMALLEST TO LARGEST) OF THE * SMALLEST EIGENVALUE TO BE RETURNED. IL MUST BE AT LEAST 1. * NOT REFERENCED IF RANGE='A' OR 'V'. * * IU (INPUT) INTEGER * IF RANGE='I', THE INDEX (FROM SMALLEST TO LARGEST) OF THE * LARGEST EIGENVALUE TO BE RETURNED. IU MUST BE AT LEAST IL * AND NO GREATER THAN N. NOT REFERENCED IF RANGE='A' OR 'V'. * * ABSTOL (INPUT) DOUBLE PRECISION * THE ABSOLUTE TOLERANCE FOR THE EIGENVALUES. AN EIGENVALUE * (OR CLUSTER) IS CONSIDERED TO BE LOCATED IF IT HAS BEEN * DETERMINED TO LIE IN AN INTERVAL WHOSE WIDTH IS ABSTOL OR * LESS. IF ABSTOL IS LESS THAN OR EQUAL TO ZERO, THEN ULP*|T| * WILL BE USED, WHERE |T| MEANS THE 1-NORM OF T. * * D (INPUT) DOUBLE PRECISION ARRAY, DIMENSION (N) * THE N DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX T. TO * AVOID OVERFLOW, THE MATRIX MUST BE SCALED SO THAT ITS LARGEST * ENTRY IS NO GREATER THAN OVERFLOW**(1/2) * UNDERFLOW**(1/4) * IN ABSOLUTE VALUE, AND FOR GREATEST ACCURACY, IT SHOULD NOT * BE MUCH SMALLER THAN THAT. * * E (INPUT) DOUBLE PRECISION ARRAY, DIMENSION (N-1) * THE (N-1) OFF-DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX T. * TO AVOID OVERFLOW, THE MATRIX MUST BE SCALED SO THAT ITS * LARGEST ENTRY IS NO GREATER THAN OVERFLOW**(1/2) * * UNDERFLOW**(1/4) IN ABSOLUTE VALUE, AND FOR GREATEST * ACCURACY, IT SHOULD NOT BE MUCH SMALLER THAN THAT. * * M (OUTPUT) INTEGER * THE ACTUAL NUMBER OF EIGENVALUES FOUND. 0 <= M <= N. * (SEE ALSO THE DESCRIPTION OF INFO=2,3.) * * NSPLIT (OUTPUT) INTEGER * THE NUMBER OF DIAGONAL BLOCKS IN THE MATRIX T. * 1 <= NSPLIT <= N. * * W (OUTPUT) DOUBLE PRECISION ARRAY, DIMENSION (N) * ON EXIT, THE FIRST M ELEMENTS OF W WILL CONTAIN THE * EIGENVALUES. (DSTEBZ MAY USE THE REMAINING N-M ELEMENTS AS * WORKSPACE.) * * IBLOCK (OUTPUT) INTEGER ARRAY, DIMENSION (N) * AT EACH ROW/COLUMN J WHERE E(J) IS ZERO OR SMALL, THE * MATRIX T IS CONSIDERED TO SPLIT INTO A BLOCK DIAGONAL * MATRIX. ON EXIT, IBLOCK(I) SPECIFIES WHICH BLOCK (FROM 1 TO * THE NUMBER OF BLOCKS) THE EIGENVALUE W(I) BELONGS TO. * (DSTEBZ MAY USE THE REMAINING N-M ELEMENTS AS WORKSPACE.) * * ISPLIT (OUTPUT) INTEGER ARRAY, DIMENSION (N) * THE SPLITTING POINTS, AT WHICH T BREAKS UP INTO SUBMATRICES. * THE FIRST SUBMATRIX CONSISTS OF ROWS/COLUMNS 1 TO ISPLIT(1), * THE SECOND OF ROWS/COLUMNS ISPLIT(1)+1 THROUGH ISPLIT(2), * ETC., AND THE NSPLIT-TH CONSISTS OF ROWS/COLUMNS * ISPLIT(NSPLIT-1)+1 THROUGH ISPLIT(NSPLIT)=N. * (ONLY THE FIRST NSPLIT ELEMENTS WILL ACTUALLY BE USED, BUT * SINCE THE USER CANNOT KNOW A PRIORI WHAT VALUE NSPLIT WILL * HAVE, N WORDS MUST BE RESERVED FOR ISPLIT.) * * WORK (WORKSPACE) DOUBLE PRECISION ARRAY, DIMENSION (4*N) * * IWORK (WORKSPACE) INTEGER ARRAY, DIMENSION (3*N) * * INFO (OUTPUT) INTEGER * = 0: SUCCESSFUL EXIT * < 0: IF INFO = -I, THE I-TH ARGUMENT HAD AN ILLEGAL VALUE * > 0: SOME OR ALL OF THE EIGENVALUES FAILED TO CONVERGE OR * WERE NOT COMPUTED: * =1 OR 3: BISECTION FAILED TO CONVERGE FOR SOME * EIGENVALUES; THESE EIGENVALUES ARE FLAGGED BY A * NEGATIVE BLOCK NUMBER. THE EFFECT IS THAT THE * EIGENVALUES MAY NOT BE AS ACCURATE AS THE * ABSOLUTE AND RELATIVE TOLERANCES. THIS IS * GENERALLY CAUSED BY UNEXPECTEDLY INACCURATE * ARITHMETIC. * =2 OR 3: RANGE='I' ONLY: NOT ALL OF THE EIGENVALUES * IL:IU WERE FOUND. * EFFECT: M < IU+1-IL * CAUSE: NON-MONOTONIC ARITHMETIC, CAUSING THE * STURM SEQUENCE TO BE NON-MONOTONIC. * CURE: RECALCULATE, USING RANGE='A', AND PICK * OUT EIGENVALUES IL:IU. IN SOME CASES, * INCREASING THE PARAMETER "FUDGE" MAY * MAKE THINGS WORK. * = 4: RANGE='I', AND THE GERSHGORIN INTERVAL * INITIALLY USED WAS TOO SMALL. NO EIGENVALUES * WERE COMPUTED. * PROBABLE CAUSE: YOUR MACHINE HAS SLOPPY * FLOATING-POINT ARITHMETIC. * CURE: INCREASE THE PARAMETER "FUDGE", * RECOMPILE, AND TRY AGAIN. * * INTERNAL PARAMETERS * =================== * * RELFAC DOUBLE PRECISION, DEFAULT = 2.0E0 * THE RELATIVE TOLERANCE. AN INTERVAL (A,B] LIES WITHIN * "RELATIVE TOLERANCE" IF B-A < RELFAC*ULP*MAX(|A|,|B|), * WHERE "ULP" IS THE MACHINE PRECISION (DISTANCE FROM 1 TO * THE NEXT LARGER FLOATING POINT NUMBER.) * * FUDGE DOUBLE PRECISION, DEFAULT = 2 * A "FUDGE FACTOR" TO WIDEN THE GERSHGORIN INTERVALS. IDEALLY, * A VALUE OF 1 SHOULD WORK, BUT ON MACHINES WITH SLOPPY * ARITHMETIC, THIS NEEDS TO BE LARGER. THE DEFAULT FOR * PUBLICLY RELEASED VERSIONS SHOULD BE LARGE ENOUGH TO HANDLE * THE WORST MACHINE AROUND. NOTE THAT THIS HAS NO EFFECT * ON ACCURACY OF THE SOLUTION. * * ===================================================================== * * .. PARAMETERS .. DOUBLE PRECISION ZERO, ONE, TWO, HALF PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, $ HALF = 1.0D0 / TWO ) DOUBLE PRECISION FUDGE, RELFAC PARAMETER ( FUDGE = 2.0D0, RELFAC = 2.0D0 ) * .. * .. LOCAL SCALARS .. LOGICAL NCNVRG, TOOFEW INTEGER IB, IBEGIN, IDISCL, IDISCU, IE, IEND, IINFO, $ IM, IN, IOFF, IORDER, IOUT, IRANGE, ITMAX, $ ITMP1, IW, IWOFF, J, JB, JDISC, JE, NB, NWL, $ NWU DOUBLE PRECISION ATOLI, BNORM, GL, GU, PIVMIN, RTOLI, SAFEMN, $ TMP1, TMP2, TNORM, ULP, WKILL, WL, WLU, WU, WUL * .. * .. LOCAL ARRAYS .. INTEGER IDUMMA( 1 ) * .. * .. EXTERNAL FUNCTIONS .. LOGICAL LSAME INTEGER ILAENV DOUBLE PRECISION DLAMCH EXTERNAL LSAME, ILAENV, DLAMCH * .. * .. EXTERNAL SUBROUTINES .. EXTERNAL DLAEBZ, XERBLA * .. * .. INTRINSIC FUNCTIONS .. INTRINSIC ABS, INT, LOG, MAX, MIN, SQRT * .. * .. EXECUTABLE STATEMENTS .. * INFO = 0 * * DECODE RANGE * IF( LSAME( RANGE, 'A' ) ) THEN IRANGE = 1 ELSE IF( LSAME( RANGE, 'V' ) ) THEN IRANGE = 2 ELSE IF( LSAME( RANGE, 'I' ) ) THEN IRANGE = 3 ELSE IRANGE = 0 END IF * * DECODE ORDER * IF( LSAME( ORDER, 'B' ) ) THEN IORDER = 2 ELSE IF( LSAME( ORDER, 'E' ) .OR. LSAME( ORDER, 'A' ) ) THEN IORDER = 1 ELSE IORDER = 0 END IF * * CHECK FOR ERRORS * IF( IRANGE.LE.0 ) THEN INFO = -1 ELSE IF( IORDER.LE.0 ) THEN INFO = -2 ELSE IF( N.LT.0 ) THEN INFO = -3 ELSE IF( IRANGE.EQ.2 .AND. VL.GE.VU ) THEN INFO = -5 ELSE IF( IRANGE.EQ.3 .AND. ( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) ) $ THEN INFO = -6 ELSE IF( IRANGE.EQ.3 .AND. ( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) ) $ THEN INFO = -7 END IF * IF( INFO.NE.0 ) THEN CALL XERBLA( 'DSTEBZ', -INFO ) RETURN END IF * * INITIALIZE ERROR FLAGS * INFO = 0 NCNVRG = .FALSE. TOOFEW = .FALSE. * * QUICK RETURN IF POSSIBLE * M = 0 IF( N.EQ.0 ) $ RETURN * * SIMPLIFICATIONS: * IF( IRANGE.EQ.3 .AND. IL.EQ.1 .AND. IU.EQ.N ) $ IRANGE = 1 * * GET MACHINE CONSTANTS * NB IS THE MINIMUM VECTOR LENGTH FOR VECTOR BISECTION, OR 0 * IF ONLY SCALAR IS TO BE DONE. * SAFEMN = DLAMCH( 'S' ) ULP = DLAMCH( 'P' ) RTOLI = ULP*RELFAC NB = ILAENV( 1, 'DSTEBZ', ' ', N, -1, -1, -1 ) IF( NB.LE.1 ) $ NB = 0 * * SPECIAL CASE WHEN N=1 * IF( N.EQ.1 ) THEN NSPLIT = 1 ISPLIT( 1 ) = 1 IF( IRANGE.EQ.2 .AND. ( VL.GE.D( 1 ) .OR. VU.LT.D( 1 ) ) ) THEN M = 0 ELSE W( 1 ) = D( 1 ) IBLOCK( 1 ) = 1 M = 1 END IF RETURN END IF * * COMPUTE SPLITTING POINTS * NSPLIT = 1 WORK( N ) = ZERO PIVMIN = ONE * DO 10 J = 2, N TMP1 = E( J-1 )**2 IF( ABS( D( J )*D( J-1 ) )*ULP**2+SAFEMN.GT.TMP1 ) THEN ISPLIT( NSPLIT ) = J - 1 NSPLIT = NSPLIT + 1 WORK( J-1 ) = ZERO ELSE WORK( J-1 ) = TMP1 PIVMIN = MAX( PIVMIN, TMP1 ) END IF 10 CONTINUE ISPLIT( NSPLIT ) = N PIVMIN = PIVMIN*SAFEMN * * COMPUTE INTERVAL AND ATOLI * IF( IRANGE.EQ.3 ) THEN * * RANGE='I': COMPUTE THE INTERVAL CONTAINING EIGENVALUES * IL THROUGH IU. * * COMPUTE GERSHGORIN INTERVAL FOR ENTIRE (SPLIT) MATRIX * AND USE IT AS THE INITIAL INTERVAL * GU = D( 1 ) GL = D( 1 ) TMP1 = ZERO * DO 20 J = 1, N - 1 TMP2 = SQRT( WORK( J ) ) GU = MAX( GU, D( J )+TMP1+TMP2 ) GL = MIN( GL, D( J )-TMP1-TMP2 ) TMP1 = TMP2 20 CONTINUE * GU = MAX( GU, D( N )+TMP1 ) GL = MIN( GL, D( N )-TMP1 ) TNORM = MAX( ABS( GL ), ABS( GU ) ) GL = GL - FUDGE*TNORM*ULP*N - FUDGE*TWO*PIVMIN GU = GU + FUDGE*TNORM*ULP*N + FUDGE*PIVMIN * * COMPUTE ITERATION PARAMETERS * ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) / $ LOG( TWO ) ) + 2 IF( ABSTOL.LE.ZERO ) THEN ATOLI = ULP*TNORM ELSE ATOLI = ABSTOL END IF * WORK( N+1 ) = GL WORK( N+2 ) = GL WORK( N+3 ) = GU WORK( N+4 ) = GU WORK( N+5 ) = GL WORK( N+6 ) = GU IWORK( 1 ) = -1 IWORK( 2 ) = -1 IWORK( 3 ) = N + 1 IWORK( 4 ) = N + 1 IWORK( 5 ) = IL - 1 IWORK( 6 ) = IU * CALL DLAEBZ( 3, ITMAX, N, 2, 2, NB, ATOLI, RTOLI, PIVMIN, D, E, $ WORK, IWORK( 5 ), WORK( N+1 ), WORK( N+5 ), IOUT, $ IWORK, W, IBLOCK, IINFO ) * IF( IWORK( 6 ).EQ.IU ) THEN WL = WORK( N+1 ) WLU = WORK( N+3 ) NWL = IWORK( 1 ) WU = WORK( N+4 ) WUL = WORK( N+2 ) NWU = IWORK( 4 ) ELSE WL = WORK( N+2 ) WLU = WORK( N+4 ) NWL = IWORK( 2 ) WU = WORK( N+3 ) WUL = WORK( N+1 ) NWU = IWORK( 3 ) END IF * IF( NWL.LT.0 .OR. NWL.GE.N .OR. NWU.LT.1 .OR. NWU.GT.N ) THEN INFO = 4 RETURN END IF ELSE * * RANGE='A' OR 'V' -- SET ATOLI * TNORM = MAX( ABS( D( 1 ) )+ABS( E( 1 ) ), $ ABS( D( N ) )+ABS( E( N-1 ) ) ) * DO 30 J = 2, N - 1 TNORM = MAX( TNORM, ABS( D( J ) )+ABS( E( J-1 ) )+ $ ABS( E( J ) ) ) 30 CONTINUE * IF( ABSTOL.LE.ZERO ) THEN ATOLI = ULP*TNORM ELSE ATOLI = ABSTOL END IF * IF( IRANGE.EQ.2 ) THEN WL = VL WU = VU END IF END IF * * FIND EIGENVALUES -- LOOP OVER BLOCKS AND RECOMPUTE NWL AND NWU. * NWL ACCUMULATES THE NUMBER OF EIGENVALUES .LE. WL, * NWU ACCUMULATES THE NUMBER OF EIGENVALUES .LE. WU * M = 0 IEND = 0 INFO = 0 NWL = 0 NWU = 0 * DO 70 JB = 1, NSPLIT IOFF = IEND IBEGIN = IOFF + 1 IEND = ISPLIT( JB ) IN = IEND - IOFF * IF( IN.EQ.1 ) THEN * * SPECIAL CASE -- IN=1 * IF( IRANGE.EQ.1 .OR. WL.GE.D( IBEGIN )-PIVMIN ) $ NWL = NWL + 1 IF( IRANGE.EQ.1 .OR. WU.GE.D( IBEGIN )-PIVMIN ) $ NWU = NWU + 1 IF( IRANGE.EQ.1 .OR. ( WL.LT.D( IBEGIN )-PIVMIN .AND. WU.GE. $ D( IBEGIN )-PIVMIN ) ) THEN M = M + 1 W( M ) = D( IBEGIN ) IBLOCK( M ) = JB END IF ELSE * * GENERAL CASE -- IN > 1 * * COMPUTE GERSHGORIN INTERVAL * AND USE IT AS THE INITIAL INTERVAL * GU = D( IBEGIN ) GL = D( IBEGIN ) TMP1 = ZERO * DO 40 J = IBEGIN, IEND - 1 TMP2 = ABS( E( J ) ) GU = MAX( GU, D( J )+TMP1+TMP2 ) GL = MIN( GL, D( J )-TMP1-TMP2 ) TMP1 = TMP2 40 CONTINUE * GU = MAX( GU, D( IEND )+TMP1 ) GL = MIN( GL, D( IEND )-TMP1 ) BNORM = MAX( ABS( GL ), ABS( GU ) ) GL = GL - FUDGE*BNORM*ULP*IN - FUDGE*PIVMIN GU = GU + FUDGE*BNORM*ULP*IN + FUDGE*PIVMIN * IF( IRANGE.GT.1 ) THEN GL = MAX( GL, WL ) GU = MIN( GU, WU ) IF( GL.GE.GU ) $ GO TO 70 END IF * * SET UP INITIAL INTERVAL * WORK( N+1 ) = GL WORK( N+IN+1 ) = GU CALL DLAEBZ( 1, 0, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN, $ D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ), $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IM, $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO ) * NWL = NWL + IWORK( 1 ) NWU = NWU + IWORK( IN+1 ) IWOFF = M - IWORK( 1 ) * * COMPUTE EIGENVALUES * ITMAX = INT( ( LOG( GU-GL+PIVMIN )-LOG( PIVMIN ) ) / $ LOG( TWO ) ) + 2 CALL DLAEBZ( 2, ITMAX, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN, $ D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ), $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IOUT, $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO ) * * COPY EIGENVALUES INTO W AND IBLOCK * USE -JB FOR BLOCK NUMBER FOR UNCONVERGED EIGENVALUES. * DO 60 J = 1, IOUT TMP1 = HALF*( WORK( J+N )+WORK( J+IN+N ) ) * * FLAG NON-CONVERGENCE. * IF( J.GT.IOUT-IINFO ) THEN NCNVRG = .TRUE. IB = -JB ELSE IB = JB END IF DO 50 JE = IWORK( J ) + 1 + IWOFF, $ IWORK( J+IN ) + IWOFF W( JE ) = TMP1 IBLOCK( JE ) = IB 50 CONTINUE 60 CONTINUE * M = M + IM END IF 70 CONTINUE * * IF RANGE='I', THEN (WL,WU) CONTAINS EIGENVALUES NWL+1,...,NWU * IF NWL+1 < IL OR NWU > IU, DISCARD EXTRA EIGENVALUES. * IF( IRANGE.EQ.3 ) THEN IM = 0 IDISCL = IL - 1 - NWL IDISCU = NWU - IU * IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN DO 80 JE = 1, M IF( W( JE ).LE.WLU .AND. IDISCL.GT.0 ) THEN IDISCL = IDISCL - 1 ELSE IF( W( JE ).GE.WUL .AND. IDISCU.GT.0 ) THEN IDISCU = IDISCU - 1 ELSE IM = IM + 1 W( IM ) = W( JE ) IBLOCK( IM ) = IBLOCK( JE ) END IF 80 CONTINUE M = IM END IF IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN * * CODE TO DEAL WITH EFFECTS OF BAD ARITHMETIC: * SOME LOW EIGENVALUES TO BE DISCARDED ARE NOT IN (WL,WLU], * OR HIGH EIGENVALUES TO BE DISCARDED ARE NOT IN (WUL,WU] * SO JUST KILL OFF THE SMALLEST IDISCL/LARGEST IDISCU * EIGENVALUES, BY SIMPLY FINDING THE SMALLEST/LARGEST * EIGENVALUE(S). * * (IF N(W) IS MONOTONE NON-DECREASING, THIS SHOULD NEVER * HAPPEN.) * IF( IDISCL.GT.0 ) THEN WKILL = WU DO 100 JDISC = 1, IDISCL IW = 0 DO 90 JE = 1, M IF( IBLOCK( JE ).NE.0 .AND. $ ( W( JE ).LT.WKILL .OR. IW.EQ.0 ) ) THEN IW = JE WKILL = W( JE ) END IF 90 CONTINUE IBLOCK( IW ) = 0 100 CONTINUE END IF IF( IDISCU.GT.0 ) THEN * WKILL = WL DO 120 JDISC = 1, IDISCU IW = 0 DO 110 JE = 1, M IF( IBLOCK( JE ).NE.0 .AND. $ ( W( JE ).GT.WKILL .OR. IW.EQ.0 ) ) THEN IW = JE WKILL = W( JE ) END IF 110 CONTINUE IBLOCK( IW ) = 0 120 CONTINUE END IF IM = 0 DO 130 JE = 1, M IF( IBLOCK( JE ).NE.0 ) THEN IM = IM + 1 W( IM ) = W( JE ) IBLOCK( IM ) = IBLOCK( JE ) END IF 130 CONTINUE M = IM END IF IF( IDISCL.LT.0 .OR. IDISCU.LT.0 ) THEN TOOFEW = .TRUE. END IF END IF * * IF ORDER='B', DO NOTHING -- THE EIGENVALUES ARE ALREADY SORTED * BY BLOCK. * IF ORDER='E' OR 'A', SORT THE EIGENVALUES FROM SMALLEST TO LARGEST * IF( IORDER.EQ.1 .AND. NSPLIT.GT.1 ) THEN DO 150 JE = 1, M - 1 IE = 0 TMP1 = W( JE ) DO 140 J = JE + 1, M IF( W( J ).LT.TMP1 ) THEN IE = J TMP1 = W( J ) END IF 140 CONTINUE * IF( IE.NE.0 ) THEN ITMP1 = IBLOCK( IE ) W( IE ) = W( JE ) IBLOCK( IE ) = IBLOCK( JE ) W( JE ) = TMP1 IBLOCK( JE ) = ITMP1 END IF 150 CONTINUE END IF * INFO = 0 IF( NCNVRG ) $ INFO = INFO + 1 IF( TOOFEW ) $ INFO = INFO + 2 RETURN * * END OF DSTEBZ * END CUT HERE............ CAT > DLAEBZ.F <<'CUT HERE............' SUBROUTINE DLAEBZ( IJOB, NITMAX, N, MMAX, MINP, NBMIN, ABSTOL, $ RELTOL, PIVMIN, D, E, E2, NVAL, AB, C, MOUT, $ NAB, WORK, IWORK, INFO ) * * -- LAPACK AUXILIARY ROUTINE (VERSION 1.1) -- * UNIV. OF TENNESSEE, UNIV. OF CALIFORNIA BERKELEY, NAG LTD., * COURANT INSTITUTE, ARGONNE NATIONAL LAB, AND RICE UNIVERSITY * MARCH 31, 1993 * * .. SCALAR ARGUMENTS .. INTEGER IJOB, INFO, MINP, MMAX, MOUT, N, NBMIN, NITMAX DOUBLE PRECISION ABSTOL, PIVMIN, RELTOL * .. * .. ARRAY ARGUMENTS .. INTEGER IWORK( * ), NAB( MMAX, * ), NVAL( * ) DOUBLE PRECISION AB( MMAX, * ), C( * ), D( * ), E( * ), E2( * ), $ WORK( * ) * .. * * PURPOSE * ======= * * DLAEBZ CONTAINS THE ITERATION LOOPS WHICH COMPUTE AND USE THE * FUNCTION N(W), WHICH IS THE COUNT OF EIGENVALUES OF A SYMMETRIC * TRIDIAGONAL MATRIX T LESS THAN OR EQUAL TO ITS ARGUMENT W. IT * PERFORMS A CHOICE OF TWO TYPES OF LOOPS: * * IJOB=1, FOLLOWED BY * IJOB=2: IT TAKES AS INPUT A LIST OF INTERVALS AND RETURNS A LIST OF * SUFFICIENTLY SMALL INTERVALS WHOSE UNION CONTAINS THE SAME * EIGENVALUES AS THE UNION OF THE ORIGINAL INTERVALS. * THE INPUT INTERVALS ARE (AB(J,1),AB(J,2)], J=1,...,MINP. * THE OUTPUT INTERVAL (AB(J,1),AB(J,2)] WILL CONTAIN * EIGENVALUES NAB(J,1)+1,...,NAB(J,2), WHERE 1 <= J <= MOUT. * * IJOB=3: IT PERFORMS A BINARY SEARCH IN EACH INPUT INTERVAL * (AB(J,1),AB(J,2)] FOR A POINT W(J) SUCH THAT * N(W(J))=NVAL(J), AND USES C(J) AS THE STARTING POINT OF * THE SEARCH. IF SUCH A W(J) IS FOUND, THEN ON OUTPUT * AB(J,1)=AB(J,2)=W. IF NO SUCH W(J) IS FOUND, THEN ON OUTPUT * (AB(J,1),AB(J,2)] WILL BE A SMALL INTERVAL CONTAINING THE * POINT WHERE N(W) JUMPS THROUGH NVAL(J), UNLESS THAT POINT * LIES OUTSIDE THE INITIAL INTERVAL. * * NOTE THAT THE INTERVALS ARE IN ALL CASES HALF-OPEN INTERVALS, * I.E., OF THE FORM (A,B] , WHICH INCLUDES B BUT NOT A . * * SEE W. KAHAN "ACCURATE EIGENVALUES OF A SYMMETRIC TRIDIAGONAL * MATRIX", REPORT CS41, COMPUTER SCIENCE DEPT., STANFORD * UNIVERSITY, JULY 21, 1966 * * NOTE: THE ARGUMENTS ARE, IN GENERAL, *NOT* CHECKED FOR UNREASONABLE * VALUES. * * ARGUMENTS * ========= * * IJOB (INPUT) INTEGER * SPECIFIES WHAT IS TO BE DONE: * = 1: COMPUTE NAB FOR THE INITIAL INTERVALS. * = 2: PERFORM BISECTION ITERATION TO FIND EIGENVALUES OF T. * = 3: PERFORM BISECTION ITERATION TO INVERT N(W), I.E., * TO FIND A POINT WHICH HAS A SPECIFIED NUMBER OF * EIGENVALUES OF T TO ITS LEFT. * OTHER VALUES WILL CAUSE DLAEBZ TO RETURN WITH INFO=-1. * * NITMAX (INPUT) INTEGER * THE MAXIMUM NUMBER OF "LEVELS" OF BISECTION TO BE * PERFORMED, I.E., AN INTERVAL OF WIDTH W WILL NOT BE MADE * SMALLER THAN 2^(-NITMAX) * W. IF NOT ALL INTERVALS * HAVE CONVERGED AFTER NITMAX ITERATIONS, THEN INFO IS SET * TO THE NUMBER OF NON-CONVERGED INTERVALS. * * N (INPUT) INTEGER * THE DIMENSION N OF THE TRIDIAGONAL MATRIX T. IT MUST BE AT * LEAST 1. * * MMAX (INPUT) INTEGER * THE MAXIMUM NUMBER OF INTERVALS. IF MORE THAN MMAX INTERVALS * ARE GENERATED, THEN DLAEBZ WILL QUIT WITH INFO=MMAX+1. * * MINP (INPUT) INTEGER * THE INITIAL NUMBER OF INTERVALS. IT MAY NOT BE GREATER THAN * MMAX. * * NBMIN (INPUT) INTEGER * THE SMALLEST NUMBER OF INTERVALS THAT SHOULD BE PROCESSED * USING A VECTOR LOOP. IF ZERO, THEN ONLY THE SCALAR LOOP * WILL BE USED. * * ABSTOL (INPUT) DOUBLE PRECISION * THE MINIMUM (ABSOLUTE) WIDTH OF AN INTERVAL. WHEN AN * INTERVAL IS NARROWER THAN ABSTOL, OR THAN RELTOL TIMES THE * LARGER (IN MAGNITUDE) ENDPOINT, THEN IT IS CONSIDERED TO BE * SUFFICIENTLY SMALL, I.E., CONVERGED. THIS MUST BE AT LEAST * ZERO. * * RELTOL (INPUT) DOUBLE PRECISION * THE MINIMUM RELATIVE WIDTH OF AN INTERVAL. WHEN AN INTERVAL * IS NARROWER THAN ABSTOL, OR THAN RELTOL TIMES THE LARGER (IN * MAGNITUDE) ENDPOINT, THEN IT IS CONSIDERED TO BE * SUFFICIENTLY SMALL, I.E., CONVERGED. NOTE: THIS SHOULD * ALWAYS BE AT LEAST RADIX*MACHINE EPSILON. * * PIVMIN (INPUT) DOUBLE PRECISION * THE MINIMUM ABSOLUTE VALUE OF A "PIVOT" IN THE STURM * SEQUENCE LOOP. THIS *MUST* BE AT LEAST MAX |E(J)**2| * * SAFE_MIN AND AT LEAST SAFE_MIN, WHERE SAFE_MIN IS AT LEAST * THE SMALLEST NUMBER THAT CAN DIVIDE ONE WITHOUT OVERFLOW. * * D (INPUT) DOUBLE PRECISION ARRAY, DIMENSION (N) * THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX T. TO AVOID * UNDERFLOW, THE MATRIX SHOULD BE SCALED SO THAT ITS LARGEST * ENTRY IS NO GREATER THAN OVERFLOW**(1/2) * UNDERFLOW**(1/4) * IN ABSOLUTE VALUE. TO ASSURE THE MOST ACCURATE COMPUTATION * OF SMALL EIGENVALUES, THE MATRIX SHOULD BE SCALED TO BE * NOT MUCH SMALLER THAN THAT, EITHER. * * E (INPUT) DOUBLE PRECISION ARRAY, DIMENSION (N) * THE OFFDIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX T IN * POSITIONS 1 THROUGH N-1. E(N) IS ARBITRARY. * TO AVOID UNDERFLOW, THE * MATRIX SHOULD BE SCALED SO THAT ITS LARGEST ENTRY IS NO * GREATER THAN OVERFLOW**(1/2) * UNDERFLOW**(1/4) IN ABSOLUTE * VALUE. TO ASSURE THE MOST ACCURATE COMPUTATION OF SMALL * EIGENVALUES, THE MATRIX SHOULD BE SCALED TO BE NOT MUCH * SMALLER THAN THAT, EITHER. * * E2 (INPUT) DOUBLE PRECISION ARRAY, DIMENSION (N) * THE SQUARES OF THE OFFDIAGONAL ELEMENTS OF THE TRIDIAGONAL * MATRIX T. E2(N) IS IGNORED. * * NVAL (INPUT/OUTPUT) INTEGER ARRAY, DIMENSION (MINP) * IF IJOB=1 OR 2, NOT REFERENCED. * IF IJOB=3, THE DESIRED VALUES OF N(W). THE ELEMENTS OF NVAL * WILL BE REORDERED TO CORRESPOND WITH THE INTERVALS IN AB. * THUS, NVAL(J) ON OUTPUT WILL NOT, IN GENERAL BE THE SAME AS * NVAL(J) ON INPUT, BUT IT WILL CORRESPOND WITH THE INTERVAL * (AB(J,1),AB(J,2)] ON OUTPUT. * * AB (INPUT/OUTPUT) DOUBLE PRECISION ARRAY, DIMENSION (MMAX,2) * THE ENDPOINTS OF THE INTERVALS. AB(J,1) IS A(J), THE LEFT * ENDPOINT OF THE J-TH INTERVAL, AND AB(J,2) IS B(J), THE * RIGHT ENDPOINT OF THE J-TH INTERVAL. THE INPUT INTERVALS * WILL, IN GENERAL, BE MODIFIED, SPLIT, AND REORDERED BY THE * CALCULATION. * * C (INPUT/WORKSPACE) DOUBLE PRECISION ARRAY, DIMENSION (MMAX) * IF IJOB=1, IGNORED. * IF IJOB=2, WORKSPACE. * IF IJOB=3, THEN ON INPUT C(J) SHOULD BE INITIALIZED TO THE * FIRST SEARCH POINT IN THE BINARY SEARCH. * * MOUT (OUTPUT) INTEGER * IF IJOB=1, THE NUMBER OF EIGENVALUES IN THE INTERVALS. * IF IJOB=2 OR 3, THE NUMBER OF INTERVALS OUTPUT. * IF IJOB=3, MOUT WILL EQUAL MINP. * * NAB (INPUT/OUTPUT) INTEGER ARRAY, DIMENSION (MMAX,2) * IF IJOB=1, THEN ON OUTPUT NAB(I,J) WILL BE SET TO N(AB(I,J)). * IF IJOB=2, THEN ON INPUT, NAB(I,J) SHOULD BE SET. IT MUST * SATISFY THE CONDITION: * N(AB(I,1)) <= NAB(I,1) <= NAB(I,2) <= N(AB(I,2)), * WHICH MEANS THAT IN INTERVAL I ONLY EIGENVALUES * NAB(I,1)+1,...,NAB(I,2) WILL BE CONSIDERED. USUALLY, * NAB(I,J)=N(AB(I,J)), FROM A PREVIOUS CALL TO DLAEBZ WITH * IJOB=1. * ON OUTPUT, NAB(I,J) WILL CONTAIN * MAX(NA(K),MIN(NB(K),N(AB(I,J)))), WHERE K IS THE INDEX OF * THE INPUT INTERVAL THAT THE OUTPUT INTERVAL * (AB(J,1),AB(J,2)] CAME FROM, AND NA(K) AND NB(K) ARE THE * THE INPUT VALUES OF NAB(K,1) AND NAB(K,2). * IF IJOB=3, THEN ON OUTPUT, NAB(I,J) CONTAINS N(AB(I,J)), * UNLESS N(W) > NVAL(I) FOR ALL SEARCH POINTS W , IN WHICH * CASE NAB(I,1) WILL NOT BE MODIFIED, I.E., THE OUTPUT * VALUE WILL BE THE SAME AS THE INPUT VALUE (MODULO * REORDERINGS -- SEE NVAL AND AB), OR UNLESS N(W) < NVAL(I) * FOR ALL SEARCH POINTS W , IN WHICH CASE NAB(I,2) WILL * NOT BE MODIFIED. NORMALLY, NAB SHOULD BE SET TO SOME * DISTINCTIVE VALUE(S) BEFORE DLAEBZ IS CALLED. * * WORK (WORKSPACE) DOUBLE PRECISION ARRAY, DIMENSION (MMAX) * WORKSPACE. * * IWORK (WORKSPACE) INTEGER ARRAY, DIMENSION (MMAX) * WORKSPACE. * * INFO (OUTPUT) INTEGER * = 0: ALL INTERVALS CONVERGED. * = 1--MMAX: THE LAST INFO INTERVALS DID NOT CONVERGE. * = MMAX+1: MORE THAN MMAX INTERVALS WERE GENERATED. * * FURTHER DETAILS * =============== * * THIS ROUTINE IS INTENDED TO BE CALLED ONLY BY OTHER LAPACK * ROUTINES, THUS THE INTERFACE IS LESS USER-FRIENDLY. IT IS INTENDED * FOR TWO PURPOSES: * * (A) FINDING EIGENVALUES. IN THIS CASE, DLAEBZ SHOULD HAVE ONE OR * MORE INITIAL INTERVALS SET UP IN AB, AND DLAEBZ SHOULD BE CALLED * WITH IJOB=1. THIS SETS UP NAB, AND ALSO COUNTS THE EIGENVALUES. * INTERVALS WITH NO EIGENVALUES WOULD USUALLY BE THROWN OUT AT * THIS POINT. ALSO, IF NOT ALL THE EIGENVALUES IN AN INTERVAL I * ARE DESIRED, NAB(I,1) CAN BE INCREASED OR NAB(I,2) DECREASED. * FOR EXAMPLE, SET NAB(I,1)=NAB(I,2)-1 TO GET THE LARGEST * EIGENVALUE. DLAEBZ IS THEN CALLED WITH IJOB=2 AND MMAX * NO SMALLER THAN THE VALUE OF MOUT RETURNED BY THE CALL WITH * IJOB=1. AFTER THIS (IJOB=2) CALL, EIGENVALUES NAB(I,1)+1 * THROUGH NAB(I,2) ARE APPROXIMATELY AB(I,1) (OR AB(I,2)) TO THE * TOLERANCE SPECIFIED BY ABSTOL AND RELTOL. * * (B) FINDING AN INTERVAL (A',B'] CONTAINING EIGENVALUES W(F),...,W(L). * IN THIS CASE, START WITH A GERSHGORIN INTERVAL (A,B). SET UP * AB TO CONTAIN 2 SEARCH INTERVALS, BOTH INITIALLY (A,B). ONE * NVAL ENTRY SHOULD CONTAIN F-1 AND THE OTHER SHOULD CONTAIN L * , WHILE C SHOULD CONTAIN A AND B, RESP. NAB(I,1) SHOULD BE -1 * AND NAB(I,2) SHOULD BE N+1, TO FLAG AN ERROR IF THE DESIRED * INTERVAL DOES NOT LIE IN (A,B). DLAEBZ IS THEN CALLED WITH * IJOB=3. ON EXIT, IF W(F-1) < W(F), THEN ONE OF THE INTERVALS -- * J -- WILL HAVE AB(J,1)=AB(J,2) AND NAB(J,1)=NAB(J,2)=F-1, WHILE * IF, TO THE SPECIFIED TOLERANCE, W(F-K)=...=W(F+R), K > 0 AND R * >= 0, THEN THE INTERVAL WILL HAVE N(AB(J,1))=NAB(J,1)=F-K AND * N(AB(J,2))=NAB(J,2)=F+R. THE CASES W(L) < W(L+1) AND * W(L-R)=...=W(L+K) ARE HANDLED SIMILARLY. * * ===================================================================== * * .. PARAMETERS .. DOUBLE PRECISION ZERO, TWO, HALF PARAMETER ( ZERO = 0.0D0, TWO = 2.0D0, $ HALF = 1.0D0 / TWO ) * .. * .. LOCAL SCALARS .. INTEGER ITMP1, ITMP2, J, JI, JIT, JP, KF, KFNEW, KL, $ KLNEW DOUBLE PRECISION TMP1, TMP2 * .. * .. INTRINSIC FUNCTIONS .. INTRINSIC ABS, MAX, MIN * .. * .. EXECUTABLE STATEMENTS .. * * CHECK FOR ERRORS * INFO = 0 IF( IJOB.LT.1 .OR. IJOB.GT.3 ) THEN INFO = -1 RETURN END IF * * INITIALIZE NAB * IF( IJOB.EQ.1 ) THEN * * COMPUTE THE NUMBER OF EIGENVALUES IN THE INITIAL INTERVALS. * MOUT = 0 DO 30 JI = 1, MINP DO 20 JP = 1, 2 TMP1 = D( 1 ) - AB( JI, JP ) IF( ABS( TMP1 ).LT.PIVMIN ) $ TMP1 = -PIVMIN NAB( JI, JP ) = 0 IF( TMP1.LE.ZERO ) $ NAB( JI, JP ) = 1 * DO 10 J = 2, N TMP1 = D( J ) - E2( J-1 ) / TMP1 - AB( JI, JP ) IF( ABS( TMP1 ).LT.PIVMIN ) $ TMP1 = -PIVMIN IF( TMP1.LE.ZERO ) $ NAB( JI, JP ) = NAB( JI, JP ) + 1 10 CONTINUE 20 CONTINUE MOUT = MOUT + NAB( JI, 2 ) - NAB( JI, 1 ) 30 CONTINUE RETURN END IF * * INITIALIZE FOR LOOP * * KF AND KL HAVE THE FOLLOWING MEANING: * INTERVALS 1,...,KF-1 HAVE CONVERGED. * INTERVALS KF,...,KL STILL NEED TO BE REFINED. * KF = 1 KL = MINP * * IF IJOB=2, INITIALIZE C. * IF IJOB=3, USE THE USER-SUPPLIED STARTING POINT. * IF( IJOB.EQ.2 ) THEN DO 40 JI = 1, MINP C( JI ) = HALF*( AB( JI, 1 )+AB( JI, 2 ) ) 40 CONTINUE END IF * * ITERATION LOOP * DO 130 JIT = 1, NITMAX * * LOOP OVER INTERVALS * IF( KL-KF+1.GE.NBMIN .AND. NBMIN.GT.0 ) THEN * * BEGIN OF PARALLEL VERSION OF THE LOOP * DO 60 JI = KF, KL * * COMPUTE N(C), THE NUMBER OF EIGENVALUES LESS THAN C * WORK( JI ) = D( 1 ) - C( JI ) IWORK( JI ) = 0 IF( WORK( JI ).LE.PIVMIN ) THEN IWORK( JI ) = 1 WORK( JI ) = MIN( WORK( JI ), -PIVMIN ) END IF * DO 50 J = 2, N WORK( JI ) = D( J ) - E2( J-1 ) / WORK( JI ) - C( JI ) IF( WORK( JI ).LE.PIVMIN ) THEN IWORK( JI ) = IWORK( JI ) + 1 WORK( JI ) = MIN( WORK( JI ), -PIVMIN ) END IF 50 CONTINUE 60 CONTINUE * IF( IJOB.LE.2 ) THEN * * IJOB=2: CHOOSE ALL INTERVALS CONTAINING EIGENVALUES. * KLNEW = KL DO 70 JI = KF, KL * * INSURE THAT N(W) IS MONOTONE * IWORK( JI ) = MIN( NAB( JI, 2 ), $ MAX( NAB( JI, 1 ), IWORK( JI ) ) ) * * UPDATE THE QUEUE -- ADD INTERVALS IF BOTH HALVES * CONTAIN EIGENVALUES. * IF( IWORK( JI ).EQ.NAB( JI, 2 ) ) THEN * * NO EIGENVALUE IN THE UPPER INTERVAL: * JUST USE THE LOWER INTERVAL. * AB( JI, 2 ) = C( JI ) * ELSE IF( IWORK( JI ).EQ.NAB( JI, 1 ) ) THEN * * NO EIGENVALUE IN THE LOWER INTERVAL: * JUST USE THE UPPER INTERVAL. * AB( JI, 1 ) = C( JI ) ELSE KLNEW = KLNEW + 1 IF( KLNEW.LE.MMAX ) THEN * * EIGENVALUE IN BOTH INTERVALS -- ADD UPPER TO * QUEUE. * AB( KLNEW, 2 ) = AB( JI, 2 ) NAB( KLNEW, 2 ) = NAB( JI, 2 ) AB( KLNEW, 1 ) = C( JI ) NAB( KLNEW, 1 ) = IWORK( JI ) AB( JI, 2 ) = C( JI ) NAB( JI, 2 ) = IWORK( JI ) ELSE INFO = MMAX + 1 END IF END IF 70 CONTINUE IF( INFO.NE.0 ) $ RETURN KL = KLNEW ELSE * * IJOB=3: BINARY SEARCH. KEEP ONLY THE INTERVAL CONTAINING * W S.T. N(W) = NVAL * DO 80 JI = KF, KL IF( IWORK( JI ).LE.NVAL( JI ) ) THEN AB( JI, 1 ) = C( JI ) NAB( JI, 1 ) = IWORK( JI ) END IF IF( IWORK( JI ).GE.NVAL( JI ) ) THEN AB( JI, 2 ) = C( JI ) NAB( JI, 2 ) = IWORK( JI ) END IF 80 CONTINUE END IF * ELSE * * END OF PARALLEL VERSION OF THE LOOP * * BEGIN OF SERIAL VERSION OF THE LOOP * KLNEW = KL DO 100 JI = KF, KL * * COMPUTE N(W), THE NUMBER OF EIGENVALUES LESS THAN W * TMP1 = C( JI ) TMP2 = D( 1 ) - TMP1 ITMP1 = 0 IF( TMP2.LE.PIVMIN ) THEN ITMP1 = 1 TMP2 = MIN( TMP2, -PIVMIN ) END IF * * A SERIES OF COMPILER DIRECTIVES TO DEFEAT VECTORIZATION * FOR THE NEXT LOOP * *$PL$ CMCHAR=' ' CDIR$ NEXTSCALAR C$DIR SCALAR CDIR$ NEXT SCALAR CVD$L NOVECTOR CDEC$ NOVECTOR CVD$ NOVECTOR *VDIR NOVECTOR *VOCL LOOP,SCALAR CIBM PREFER SCALAR *$PL$ CMCHAR='*' * DO 90 J = 2, N TMP2 = D( J ) - E2( J-1 ) / TMP2 - TMP1 IF( TMP2.LE.PIVMIN ) THEN ITMP1 = ITMP1 + 1 TMP2 = MIN( TMP2, -PIVMIN ) END IF 90 CONTINUE * IF( IJOB.LE.2 ) THEN * * IJOB=2: CHOOSE ALL INTERVALS CONTAINING EIGENVALUES. * * INSURE THAT N(W) IS MONOTONE * ITMP1 = MIN( NAB( JI, 2 ), $ MAX( NAB( JI, 1 ), ITMP1 ) ) * * UPDATE THE QUEUE -- ADD INTERVALS IF BOTH HALVES *