#! /bin/sh # This is a shell archive. Remove anything before this line, then unpack # it by saving it into a file and typing "sh file". To overwrite existing # files, type "sh file -c". You can also feed this as standard input via # unshar, or by typing "sh 'd1mach.f' <<'END_OF_FILE' X X double precision function d1mach( i ) X* X* -- lapack auxiliary routine -- X* argonne national lab, courant institute, and n.a.g. ltd. X* april 1, 1989 X* X* .. scalar arguments .. X integer i X* .. X* X* purpose X* ======= X* X* d1mach determines double precision machine constants by a call X* to the double precision version of machar. X* X* (see w. j. cody, "machar: a subroutine to dynamically determine X* machine parameters," toms 14, december, 1988. ) X* X* arguments X* ========= X* X* i - integer X* on entry, i is the index to one of the machine constants, X* as follows: X* X* d1mach(1) = b**(emin-1), the smallest positive magnitude X* X* d1mach(2) = b**emax*(1 - b**(-t)), the largest magnitude X* X* d1mach(3) = b**(-t), the smallest relative spacing X* X* d1mach(4) = b**(1-t), the largest relative spacing X* X* d1mach(5) = log10(b) X* X* .. local scalars .. X logical dejavu X integer ibdig, iexp, iradix, iround, machep, maxexp, X $ minexp, negeps, nguard X double precision eps, epsneg, xmax, xmin X* .. X* .. local arrays .. X double precision dmach( 5 ) X* .. X* .. external subroutines .. X external dmachr X* .. X* .. intrinsic functions .. X intrinsic dble, log10 X* .. X* .. save statement .. X save X* .. X* .. data statements .. X data dejavu / .false. / X* .. X* .. executable statements .. X* X if( .not.dejavu ) then X call dmachr( iradix, ibdig, iround, nguard, machep, negeps, X $ iexp, minexp, maxexp, eps, epsneg, xmin, xmax ) X dmach( 1 ) = xmin X dmach( 2 ) = xmax X dmach( 3 ) = eps X dmach( 4 ) = dble( iradix )*eps X dmach( 5 ) = log10( dble( iradix ) ) X dejavu = .true. X end if X* X if( i.lt.1 .or. i.gt.5 ) then X write( *, fmt = 9999 )i X 9999 format( ' d1mach - i out of bounds', i10 ) X stop X else X d1mach = dmach( i ) X end if X* X return X* X* end of d1mach X* X end END_OF_FILE if test 2192 -ne `wc -c <'d1mach.f'`; then echo shar: \"'d1mach.f'\" unpacked with wrong size! fi # end of 'd1mach.f' fi if test -f 'daxpy.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'daxpy.f'\" else echo shar: Extracting \"'daxpy.f'\" \(1556 characters\) sed "s/^X//" >'daxpy.f' <<'END_OF_FILE' X SUBROUTINE DAXPY( N, DA, DX, INCX, DY, INCY ) X* X* constant times a vector plus a vector. X* uses unrolled loops for increments equal to one. X* jack dongarra, linpack, 3/11/78. X* X* .. Scalar Arguments .. X INTEGER INCX, INCY, N X DOUBLE PRECISION DA X* .. X* .. Array Arguments .. X DOUBLE PRECISION DX( 1 ), DY( 1 ) X* .. X* .. Local Scalars .. X INTEGER I, IX, IY, M, MP1 X* .. X* .. Intrinsic Functions .. X INTRINSIC MOD X* .. X* .. Executable Statements .. X* X IF( N.LE.0 ) X $ RETURN X IF( DA.EQ.0.0D0 ) X $ RETURN X IF( INCX.EQ.1 .AND. INCY.EQ.1 ) X $ GO TO 20 X* X* code for unequal increments or equal increments X* not equal to 1 X* X IX = 1 X IY = 1 X IF( INCX.LT.0 ) X $ IX = ( -N+1 )*INCX + 1 X IF( INCY.LT.0 ) X $ IY = ( -N+1 )*INCY + 1 X DO 10 I = 1, N X DY( IY ) = DY( IY ) + DA*DX( IX ) X IX = IX + INCX X IY = IY + INCY X 10 CONTINUE X RETURN X* X* code for both increments equal to 1 X* X* X* clean-up loop X* X 20 M = MOD( N, 4 ) X IF( M.EQ.0 ) X $ GO TO 40 X DO 30 I = 1, M X DY( I ) = DY( I ) + DA*DX( I ) X 30 CONTINUE X IF( N.LT.4 ) X $ RETURN X 40 MP1 = M + 1 X DO 50 I = MP1, N, 4 X DY( I ) = DY( I ) + DA*DX( I ) X DY( I+1 ) = DY( I+1 ) + DA*DX( I+1 ) X DY( I+2 ) = DY( I+2 ) + DA*DX( I+2 ) X DY( I+3 ) = DY( I+3 ) + DA*DX( I+3 ) X 50 CONTINUE X RETURN X END END_OF_FILE if test 1556 -ne `wc -c <'daxpy.f'`; then echo shar: \"'daxpy.f'\" unpacked with wrong size! fi # end of 'daxpy.f' fi if test -f 'dchk21.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'dchk21.f'\" else echo shar: Extracting \"'dchk21.f'\" \(32890 characters\) sed "s/^X//" >'dchk21.f' <<'END_OF_FILE' X SUBROUTINE DCHK21( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, X $ NOUNIT, A, LDA, H, T1, T2, U, LDU, Z, UZ, WR1, X $ WI1, WR3, WI3, EVECTL, EVECTR, EVECTY, EVECTX, X $ WORK, NWORK, IWORK, LWORK, RESULT, INFO ) X* X* -- LAPACK test routine (preliminary version) -- X* Univ. of Tennessee, Oak Ridge National Lab, Argonne National Lab, X* Courant Institute, NAG Ltd., and Rice University X* March 26, 1990 X* X* .. Scalar Arguments .. X* X***************************** X INTEGER INFO, LDA, LDU, NOUNIT, NSIZES, NTYPES, NWORK X DOUBLE PRECISION THRESH X* .. X* X* .. Array Arguments .. X* X LOGICAL DOTYPE( * ), LWORK( * ) X INTEGER ISEED( 4 ), IWORK( * ), NN( * ) X DOUBLE PRECISION A( LDA, * ), EVECTL( LDU, * ), X $ EVECTR( LDU, * ), EVECTX( LDU, * ), X $ EVECTY( LDU, * ), H( LDA, * ), RESULT( 12 ), X $ T1( LDA, * ), T2( LDA, * ), U( LDU, * ), X $ UZ( LDU, * ), WI1( * ), WI3( * ), WORK( * ), X $ WR1( * ), WR3( * ), Z( LDU, * ) X* .. X* X*----------------------------------------------------------------------- X* X* Purpose X* ======= X* X* DCHK21 checks the nonsymmetric eigenvalue problem routines. X* X* DGEHD3 factors A as U H U' , where ' means transpose, X* H is hessenberg, and U is an orthogonal matrix. X* X* DHSEQR factors H as Z T Z' , where Z is orthogonal and X* T is "quasi-triangular", and the eigenvalue vector W. X* X* DTREVC computes the left and right eigenvector matrices X* L and R for T. X* X* DHSEIN computes the left and right eigenvector matrices X* Y and X for H, using inverse iteration. X* X* When DCHK21 is called, a number of matrix "sizes" ("n's") and a X* number of matrix "types" are specified. For each size ("n") X* and each type of matrix, one matrix will be generated and used X* to test the nonsymmetric eigenroutines. For each matrix, 12 X* tests will be performed: X* X* X* (1) | A - U H U' | / ( |A| n ulp ) X* X* (2) | I - UU' | / ( n ulp ) X* X* (3) | H - Z T Z' | / ( |H| n ulp ) X* X* (4) | I - ZZ' | / ( n ulp ) X* X* (5) | A - UZ H (UZ)' | / ( |A| n ulp ) X* X* (6) | I - UZ (UZ)' | / ( n ulp ) X* X* (7) | T(Z computed) - T(Z not computed) | / ( |T| ulp ) X* X* (8) | W(Z computed) - W(Z not computed) | / ( |W| ulp ) X* X* (9) | TR - RW | / ( |T| |R| ulp ) X* X* (10) | LT - WL | / ( |T| |L| ulp ) X* X* (11) | HX - XW | / ( |H| |X| ulp ) X* X* (12) | YH - WY | / ( |H| |Y| ulp ) X* X* The "sizes" are specified by an array NN(1:NSIZES); the value of X* each element NN(j) specifies one size. X* The "types" are specified by a logical array DOTYPE( 1:NTYPES ); X* if DOTYPE(j) is .TRUE., then matrix type "j" will be generated. X* Currently, the list of possible types is: X* X* (1) The zero matrix. X* (2) The identity matrix. X* (3) A (transposed) Jordan block, with 1's on the diagonal. X* X* (4) A diagonal matrix with evenly spaced entries X* 1, ..., ULP and random signs. X* (ULP = (first number larger than 1) - 1 ) X* (5) A diagonal matrix with geometrically spaced entries X* 1, ..., ULP and random signs. X* (6) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP X* and random signs. X* X* (7) Same as (4), but multiplied by SQRT( overflow threshold ) X* (8) Same as (4), but multiplied by SQRT( underflow threshold ) X* X* (9) A matrix of the form U' T U, where U is orthogonal and X* T has evenly spaced entries 1, ..., ULP with random signs X* on the diagonal and random O(1) entries in the upper X* triangle. X* X* (10) A matrix of the form U' T U, where U is orthogonal and X* T has geometrically spaced entries 1, ..., ULP with random X* signs on the diagonal and random O(1) entries in the upper X* triangle. X* X* (11) A matrix of the form U' T U, where U is orthogonal and X* T has "clustered" entries 1, ULP,..., ULP with random X* signs on the diagonal and random O(1) entries in the upper X* triangle. X* X* (12) A matrix of the form U' T U, where U is orthogonal and X* T has real or complex conjugate paired eigenvalues randomly X* chosen from ( ULP, 1 ) and random O(1) entries in the upper X* triangle. X* X* (13) A matrix of the form X' T X, where X has condition X* SQRT( ULP ) and T has evenly spaced entries 1, ..., ULP X* with random signs on the diagonal and random O(1) entries X* in the upper triangle. X* X* (14) A matrix of the form X' T X, where X has condition X* SQRT( ULP ) and T has geometrically spaced entries X* 1, ..., ULP with random signs on the diagonal and random X* O(1) entries in the upper triangle. X* X* (15) A matrix of the form X' T X, where X has condition X* SQRT( ULP ) and T has "clustered" entries 1, ULP,..., ULP X* with random signs on the diagonal and random O(1) entries X* in the upper triangle. X* X* (16) A matrix of the form X' T X, where X has condition X* SQRT( ULP ) and T has real or complex conjugate paired X* eigenvalues randomly chosen from ( ULP, 1 ) and random X* O(1) entries in the upper triangle. X* X* (17) Same as (16), but multiplied by SQRT( overflow threshold ) X* (18) Same as (16), but multiplied by SQRT( underflow threshold ) X* X* (19) Nonsymmetric matrix with random entries chosen from (-1,1). X* (20) Same as (19), but multiplied by SQRT( overflow threshold ) X* (21) Same as (19), but multiplied by SQRT( underflow threshold ) X* X* X* Arguments X* ========== X* X* NSIZES - INTEGER X* The number of sizes of matrices to use. If it is zero, X* DCHK21 does nothing. It must be at least zero. X* Not modified. X* X* NN - INTEGER array of dimension ( NSIZES ) X* An array containing the sizes to be used for the matrices. X* Zero values will be skipped. The values must be at least X* zero. X* Not modified. X* X* NTYPES - INTEGER X* The number of elements in DOTYPE. If it is zero, DCHK21 X* does nothing. It must be at least zero. If it is MAXTYP+1 X* and NSIZES is 1, then an additional type, MAXTYP+1 is X* defined, which is to use whatever matrix is in A. This X* is only useful if DOTYPE(1:MAXTYP) is .FALSE. and X* DOTYPE(MAXTYP+1) is .TRUE. . X* Not modified. X* X* DOTYPE - LOGICAL array of dimension ( NTYPES ) X* If DOTYPE(j) is .TRUE., then for each size in NN a X* matrix of that size and of type j will be generated. X* If NTYPES is smaller than the maximum number of types X* defined (PARAMETER MAXTYP), then types NTYPES+1 through X* MAXTYP will not be generated. If NTYPES is larger X* than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES) X* will be ignored. X* Not modified. X* X* THRESH - DOUBLE PRECISION X* A test will count as "failed" if the "error", computed as X* described above, exceeds THRESH. Note that the error X* is scaled to be O(1), so THRESH should be a reasonably X* small multiple of 1, e.g., 10 or 100. In particular, X* it should not depend on the precision (single vs. double) X* or the size of the matrix. It must be at least zero. X* Not modified. X* X* ISEED - INTEGER array of dimension ( 4 ) X* On entry ISEED specifies the seed of the random number X* generator. The array elements should be between 0 and 4095; X* if not they will be reduced mod 4096. Also, ISEED(4) must X* be odd. The random number generator uses a linear X* congruential sequence limited to small integers, and so X* should produce machine independent random numbers. The X* values of ISEED are changed on exit, and can be used in the X* next call to DCHK21 to continue the same random number X* sequence. X* Modified. X* X* NOUNIT - INTEGER X* The FORTRAN unit number for printing out error messages X* (e.g., if a routine returns INFO not equal to 0.) X* Not modified. X* X* A - DOUBLE PRECISION array of dimension ( LDA , max(NN) ) X* Used to hold the matrix whose eigenvalues are to be X* computed. On exit, A contains the last matrix actually X* used. X* Modified. X* X* LDA - INTEGER X* The leading dimension of A, H, T1, and T2. It must be at X* least 1 and at least max( NN ). X* Not modified. X* X* H - DOUBLE PRECISION array of dimension( LDA , max(NN) ) X* The upper hessenberg matrix computed by DGEHD3. On exit, X* H contains the Hessenberg form of the matrix in A. X* Modified. X* X* T1 - DOUBLE PRECISION array of dimension( LDA , max(NN) ) X* The Schur (="quasi-triangular") matrix computed by DHSEQR X* if Z is computed. On exit, T1 contains the Schur form of X* the matrix in A. X* Modified. X* X* T2 - DOUBLE PRECISION array of dimension( LDA , max(NN) ) X* The Schur matrix computed by DHSEQR when Z is not computed. X* This should be identical to T1. X* Modified. X* X* LDU - INTEGER X* The leading dimension of U, V, Z, and UZ. It must be at X* least 1 and at least max( NN ). X* Not modified. X* X* U - DOUBLE PRECISION array of dimension ( LDU, max(NN) ). X* The orthogonal matrix computed by DGEHD3. X* Modified. X* X* Z - DOUBLE PRECISION array of dimension ( LDU, max(NN) ). X* The orthogonal matrix computed by DHSEQR. X* Modified. X* X* UZ - DOUBLE PRECISION array of dimension ( LDU, max(NN) ). X* The product of U times Z. X* Modified. X* X* WR1, WI1 - DOUBLE PRECISION arrays of dimension ( max(NN) ). X* The real and imaginary parts of the eigenvalues of A, X* as computed when Z is computed. X* On exit, WR1 + WI1*i are the eigenvalues of the matrix in A. X* Modified. X* X* WR3, WI3 - DOUBLE PRECISION arrays of dimension ( max(NN) ). X* Like WR1, WI1, these arrays contain the eigenvalues of A, X* but those computed when DHSEQR only computes the X* eigenvalues, i.e., not the Schur vectors and no more of the X* Schur form than is necessary for computing the X* eigenvalues. X* Modified. X* X* EVECTL - DOUBLE PRECISION array of dimension ( LDU, max(NN) ). X* The (upper triangular) left eigenvector matrix for the X* matrix in T1. For complex conjugate pairs, the real part X* is stored in one row and the imaginary part in the next. X* Modified. X* X* EVEZTR - DOUBLE PRECISION array of dimension ( LDU, max(NN) ). X* The (upper triangular) right eigenvector matrix for the X* matrix in T1. For complex conjugate pairs, the real part X* is stored in one column and the imaginary part in the next. X* Modified. X* X* EVECTY - DOUBLE PRECISION array of dimension ( LDU, max(NN) ). X* The left eigenvector matrix for the X* matrix in H. For complex conjugate pairs, the real part X* is stored in one row and the imaginary part in the next. X* Modified. X* X* EVECTX - DOUBLE PRECISION array of dimension ( LDU, max(NN) ). X* The right eigenvector matrix for the X* matrix in H. For complex conjugate pairs, the real part X* is stored in one column and the imaginary part in the next. X* Modified. X* X* WORK - DOUBLE PRECISION array of dimension ( NWORK ) X* Workspace. X* Modified. X* X* NWORK - INTEGER X* The number of entries in WORK. This must be at least X* NN(j) * MAX( 3*NBLOCK+NSHIFT+2 , 2*NBLOCK+2*NSHIFT+2 , X* 2*NN(j) ). X* In this formula, "NBLOCK" is the blocksize and "NSHIFT" is X* the number of simultaneous shifts; both are returned by X* "ENVIR"; NBLOCK will be constrained to be between 1 and X* NN(j), while NSHIFT will be constrained to be between 2 and X* NN(j). X* Not modified. X* X* IWORK - INTEGER array of dimension ( max(NN) ) (scratch) X* Workspace. X* Modified. X* X* LWORK - LOGICAL array of dimension ( max(NN) ) (scratch) X* Workspace. Could be equivalenced to IWORK. X* Modified. X* X* RESULT - DOUBLE PRECISION array of dimension ( 12 ) (OUTPUT) X* The values computed by the twelve tests described above. X* The values are currently limited to 1/ulp, to avoid X* overflow. X* Modified. X* X* INFO - INTEGER X* If 0, then everything ran OK. X* -1: NSIZES < 0 X* -2: Some NN(j) < 0 X* -3: NTYPES < 0 X* -6: THRESH < 0 X* -9: LDA < 1 or LDA < NMAX, where NMAX is max( NN(j) ). X* -14: LDU < 1 or LDU < NMAX. X* -26: NWORK too small. X* If DLATMR, SLATMS, or SLATME returns an error code, the X* absolute value of it is returned. X* If 1, then DHSEQR could not find all the shifts. X* If 2, then the EISPACK code (for small blocks) failed. X* If >2, then 30*N iterations were not enough to find an X* eigenvalue or to decompose the problem. X* Modified. X* X* X* X*----------------------------------------------------------------------- X* X* Some Local Variables and Parameters: X* ---- ----- --------- --- ---------- X* X* ZERO, ONE Real 0 and 1. X* MAXTYP The number of types defined. X* MTEST The number of tests defined: care must be taken X* that (1) the size of RESULT, (2) the number of X* tests actually performed, and (3) MTEST agree. X* NBLOCK, NSHIFT Blocksize and number of shifts as returned by X* ENVIR. X* NMAX Largest value in NN. X* NMATS The number of matrices generated so far. X* NERRS The number of tests which have exceeded THRESH X* so far (computed by DLAFTS). X* COND, CONDS, X* IMODE Values to be passed to the matrix generators. X* ANORM Norm of A; passed to matrix generators. X* X* OVFL, UNFL Overflow and underflow thresholds. X* ULP, ULPINV Finest relative precision and its inverse. X* RTOVFL, RTUNFL, X* RTULP, RTULPI Square roots of the previous 4 values. X* X* The following four arrays decode JTYPE: X* KTYPE(j) The general type (1-10) for type "j". X* KMODE(j) The MODE value to be passed to the matrix X* generator for type "j". X* KMAGN(j) The order of magnitude ( O(1), X* O(overflow^(1/2) ), O(underflow^(1/2) ) X* KCONDS(j) Selectw whether CONDS is to be 1 or X* 1/sqrt(ulp). (0 means irrelevant.) X* X*----------------------------------------------------------------------- X* X* X* .. Parameters .. X* X DOUBLE PRECISION ZERO, ONE X PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) X INTEGER MAXTYP X PARAMETER ( MAXTYP = 21 ) X* .. X* X* .. Local Scalars .. X* X LOGICAL BADNN X INTEGER IINFO, IMODE, IN, ITYPE, J, JCOL, JSIZE, X $ JTYPE, MTYPES, N, N1, NBLOCK, NCWORK, NERRS, X $ NMATS, NMAX, NSHIFT, NTEST, NTESTT X DOUBLE PRECISION ANINV, ANORM, COND, CONDS, OVFL, RTOVFL, RTULP, X $ RTULPI, RTUNFL, TEMP1, TEMP2, ULP, ULPINV, UNFL X* .. X* X* .. Local Arrays .. X* X CHARACTER ADUMMA( 1 ) X INTEGER IDUMMA( 1 ), IOLDSD( 4 ), KCONDS( MAXTYP ), X $ KMAGN( MAXTYP ), KMODE( MAXTYP ), X $ KTYPE( MAXTYP ) X DOUBLE PRECISION DUMMA( 6 ) X* .. X* X* .. External Functions .. X* X DOUBLE PRECISION DLAMCH X EXTERNAL DLAMCH X* .. X* X* .. External Subroutines .. X* X EXTERNAL ENVIR, DGEHD3, DGEMM, DGET21, DGET22, DHSEIN, X $ DHSEQR, DLABAD, DLACPY, DLAFTS, DLASUM, DLATME, X $ DLATMR, DLATMS, DLAZRO, DTREVC, XERBLA X* .. X* X* .. Intrinsic Functions .. X* X INTRINSIC ABS, DBLE, MAX, MIN, SQRT X* .. X* X* .. Data statements .. X* X* X* X* X DATA KTYPE / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 / X DATA KMAGN / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2, X $ 3, 1, 2, 3 / X DATA KMODE / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3, X $ 1, 5, 5, 5, 4, 3, 1 / X DATA KCONDS / 3*0, 5*0, 4*1, 6*2, 3*0 / X* .. X* X* X*----------------------------------------------------------------------- X* .. Executable Statements .. X* X* X* Check for errors X* X* X* X NTESTT = 0 X INFO = 0 X* X* Important constants X* X CALL ENVIR( 'B', NBLOCK ) X CALL ENVIR( 'S', NSHIFT ) X* X BADNN = .FALSE. X NMAX = 0 X DO 10 J = 1, NSIZES X NMAX = MAX( NMAX, NN( J ) ) X IF( NN( J ).LT.0 ) X $ BADNN = .TRUE. X 10 CONTINUE X* X NBLOCK = MAX( 1, MIN( NMAX, NBLOCK ) ) X NSHIFT = MAX( 2, MIN( NMAX, NSHIFT ) ) X* X* X* Check for errors X* X* X IF( NSIZES.LT.0 ) THEN X INFO = -1 X ELSE IF( BADNN ) THEN X INFO = -2 X ELSE IF( NTYPES.LT.0 ) THEN X INFO = -3 X ELSE IF( THRESH.LT.ZERO ) THEN X INFO = -6 X ELSE IF( LDA.LE.1 .OR. LDA.LT.NMAX ) THEN X INFO = -9 X ELSE IF( LDU.LE.1 .OR. LDU.LT.NMAX ) THEN X INFO = -14 X ELSE IF( NMAX*MAX( 3*NBLOCK+NSHIFT+2, 2*NBLOCK+2*NSHIFT+2 ).GT. X $ NWORK ) THEN X INFO = -26 X END IF X* X IF( INFO.NE.0 ) THEN X CALL XERBLA( 'DCHK21', -INFO ) X RETURN X END IF X* X* Quick return if nothing to do X* X IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 ) X $ RETURN X* X* More Important constants X* X* X UNFL = DLAMCH( 'Safe minimum' ) X OVFL = DLAMCH( 'Overflow' ) X CALL DLABAD( UNFL, OVFL ) X ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' ) X ULPINV = ONE / ULP X RTUNFL = SQRT( UNFL ) X RTOVFL = SQRT( OVFL ) X RTULP = SQRT( ULP ) X RTULPI = ONE / RTULP X* X*----------------------------------------------------------------------- X* X* Loop over sizes, types X* X NERRS = 0 X NMATS = 0 X* X DO 180 JSIZE = 1, NSIZES X N = NN( JSIZE ) X N1 = MAX( 1, N ) X NCWORK = NWORK / N1 X ANINV = ONE / DBLE( N1 ) X* X* X* X IF( NSIZES.NE.1 ) THEN X MTYPES = MIN( MAXTYP, NTYPES ) X ELSE X MTYPES = MIN( MAXTYP+1, NTYPES ) X END IF X* X DO 170 JTYPE = 1, MTYPES X IF( .NOT.DOTYPE( JTYPE ) ) X $ GO TO 170 X NMATS = NMATS + 1 X NTEST = 0 X* X* Save ISEED in case of an error. X* X DO 20 J = 1, 4 X IOLDSD( J ) = ISEED( J ) X 20 CONTINUE X* X* Initialize RESULT X* X DO 30 J = 1, 12 X RESULT( J ) = ZERO X 30 CONTINUE X* X*----------------------------------------------------------------------- X* X* X* Compute "A" X* X* Control parameters: X* X* KMAGN KCONDS KMODE KTYPE X* =1 O(1) 1 clustered 1 zero X* =2 large large clustered 2 identity X* =3 small exponential Jordan X* =4 arithmetic diagonal, (w/ eigenvalues) X* =5 random log symmetric, w/ eigenvalues X* =6 random general, w/ eigenvalues X* =7 random diagonal X* =8 random symmetric X* =9 random general X* =10 random triangular X* X* X* X IF( MTYPES.GT.MAXTYP ) X $ GO TO 100 X* X ITYPE = KTYPE( JTYPE ) X IMODE = KMODE( JTYPE ) X* X* Compute norm X* X GO TO ( 40, 50, 60 )KMAGN( JTYPE ) X* X 40 CONTINUE X ANORM = ONE X GO TO 70 X* X 50 CONTINUE X ANORM = ( RTOVFL*ULP )*ANINV X GO TO 70 X* X 60 CONTINUE X ANORM = RTUNFL*N*ULPINV X GO TO 70 X* X 70 CONTINUE X* X CALL DLAZRO( LDA, N, ZERO, ZERO, A, LDA ) X IINFO = 0 X COND = ULPINV X* X* Special Matrices -- Identity & Jordan block X* X* Zero X* X IF( ITYPE.EQ.1 ) THEN X IINFO = 0 X* X* Identity X* X ELSE IF( ITYPE.EQ.2 ) THEN X* X DO 80 JCOL = 1, N X A( JCOL, JCOL ) = ANORM X 80 CONTINUE X* X* Jordan Block X* X ELSE IF( ITYPE.EQ.3 ) THEN X* X DO 90 JCOL = 1, N X A( JCOL, JCOL ) = ANORM X IF( JCOL.GT.1 ) X $ A( JCOL, JCOL-1 ) = ONE X 90 CONTINUE X* X* X* Diagonal Matrix, [Eigen]values Specified X* X ELSE IF( ITYPE.EQ.4 ) THEN X* X CALL DLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND, X $ ANORM, 0, 0, 'N', A, LDA, WORK( N+1 ), X $ IINFO ) X* X* X* Symmetric, eigenvalues specified X* X* X ELSE IF( ITYPE.EQ.5 ) THEN X* X CALL DLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND, X $ ANORM, N, N, 'N', A, LDA, WORK( N+1 ), X $ IINFO ) X* X* X* General, eigenvalues specified X* X* X ELSE IF( ITYPE.EQ.6 ) THEN X* X IF( KCONDS( JTYPE ).EQ.1 ) THEN X CONDS = ONE X ELSE IF( KCONDS( JTYPE ).EQ.2 ) THEN X CONDS = RTULPI X ELSE X CONDS = ZERO X END IF X ADUMMA( 1 ) = ' ' X CALL DLATME( N, 'S', ISEED, WORK, IMODE, COND, ONE, X $ ADUMMA, 'T', 'T', 'T', WORK( N+1 ), 4, X $ CONDS, N, N, ANORM, A, LDA, WORK( 2*N+1 ), X $ IINFO ) X* X* X* Diagonal, random eigenvalues X* X* X ELSE IF( ITYPE.EQ.7 ) THEN X* X CALL DLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE, X $ 'T', 'N', WORK( N+1 ), 1, ONE, X $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0, X $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) X* X* X* Symmetric, random eigenvalues X* X ELSE IF( ITYPE.EQ.8 ) THEN X* X CALL DLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE, X $ 'T', 'N', WORK( N+1 ), 1, ONE, X $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N, X $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) X* X* X* General, random eigenvalues X* X* X ELSE IF( ITYPE.EQ.9 ) THEN X* X CALL DLATMR( N, N, 'S', ISEED, 'N', WORK, 6, ONE, ONE, X $ 'T', 'N', WORK( N+1 ), 1, ONE, X $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N, X $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) X* X* X* Triangular, random eigenvalues X* X* X ELSE IF( ITYPE.EQ.10 ) THEN X* X CALL DLATMR( N, N, 'S', ISEED, 'N', WORK, 6, ONE, ONE, X $ 'T', 'N', WORK( N+1 ), 1, ONE, X $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, 0, X $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) X* X ELSE X IINFO = 1 X END IF X* X IF( IINFO.NE.0 ) THEN X WRITE( NOUNIT, FMT = 9999 )'Generator', IINFO, N, JTYPE, X $ IOLDSD X INFO = ABS( IINFO ) X RETURN X END IF X* X 100 CONTINUE X* X*----------------------------------------------------------------------- X* X* X* Call DGEHD3 to compute H and U, do tests. X* X* X CALL DLACPY( ' ', N, N, A, LDA, H, LDA ) X* X NTEST = 1 X CALL DGEHD3( 'Q', N, H, LDA, U, LDU, WORK, WORK( N+1 ), N1, X $ NCWORK-1, IINFO ) X IF( IINFO.NE.0 ) THEN X RESULT( 1 ) = ULPINV X WRITE( NOUNIT, FMT = 9999 )'DGEHD3', IINFO, N, JTYPE, X $ IOLDSD X INFO = ABS( IINFO ) X GO TO 160 X END IF X NTEST = 2 X* X* X CALL DGET21( 1, N, A, LDA, H, LDA, U, LDU, DUMMA, WORK, X $ RESULT( 1 ) ) X* X*----------------------------------------------------------------------- X* X* X* Call DHSEQR to compute T1, T2, and Z, do tests. X* X* X* Compute T1 and UZ X* X CALL DLACPY( ' ', N, N, H, LDA, T2, LDA ) X NTEST = 3 X RESULT( 3 ) = ULPINV X* X CALL DHSEQR( 'E', N, T2, LDA, UZ, LDU, WR3, WI3, WORK, X $ NWORK, IINFO ) X IF( IINFO.NE.0 ) THEN X WRITE( NOUNIT, FMT = 9999 )'DHSEQR(E)', IINFO, N, JTYPE, X $ IOLDSD X IF( IINFO.LE.N+2 ) THEN X INFO = ABS( IINFO ) X GO TO 160 X END IF X END IF X* X CALL DLACPY( ' ', N, N, H, LDA, T2, LDA ) X* X CALL DHSEQR( 'S', N, T2, LDA, UZ, LDU, WR1, WI1, WORK, X $ NWORK, IINFO ) X IF( IINFO.NE.0 .AND. IINFO.LE.N+2 ) THEN X WRITE( NOUNIT, FMT = 9999 )'DHSEQR(S)', IINFO, N, JTYPE, X $ IOLDSD X INFO = ABS( IINFO ) X GO TO 160 X END IF X* X CALL DLACPY( ' ', N, N, H, LDA, T1, LDA ) X CALL DLACPY( ' ', N, N, U, LDU, UZ, LDA ) X* X CALL DHSEQR( 'V', N, T1, LDA, UZ, LDU, WR1, WI1, WORK, X $ NWORK, IINFO ) X IF( IINFO.NE.0 .AND. IINFO.LE.N+2 ) THEN X WRITE( NOUNIT, FMT = 9999 )'DHSEQR(V)', IINFO, N, JTYPE, X $ IOLDSD X INFO = ABS( IINFO ) X GO TO 160 X END IF X* X* Compute Z = U' UZ X* X CALL DGEMM( 'C', 'N', N, N, N, ONE, U, LDU, UZ, LDU, ZERO, X $ Z, LDU ) X NTEST = 8 X* X* X* Do Tests 3 and 4 X* X* X CALL DGET21( 1, N, H, LDA, T1, LDA, Z, LDU, DUMMA, WORK, X $ RESULT( 3 ) ) X* X* X* X* Do Tests 5 & 6 X* X* X CALL DGET21( 1, N, A, LDA, T1, LDA, UZ, LDU, DUMMA, WORK, X $ RESULT( 5 ) ) X* X* X* Do Test 7 X* X CALL DGET21( 2, N, T2, LDA, T1, LDA, DUMMA, LDU, DUMMA, X $ WORK, RESULT( 7 ) ) X* X* X* Do Test 8 X* X TEMP1 = ZERO X TEMP2 = ZERO X* X DO 110 J = 1, N X TEMP1 = MAX( TEMP1, ABS( WR1( J ) )+ABS( WI1( J ) ), X $ ABS( WR3( J ) )+ABS( WI3( J ) ) ) X TEMP2 = MAX( TEMP2, ABS( WR1( J )-WR3( J ) )+ X $ ABS( WR1( J )-WR3( J ) ) ) X 110 CONTINUE X* X RESULT( 8 ) = TEMP2 / MAX( UNFL, ULP*MAX( TEMP1, TEMP2 ) ) X* X* X*----------------------------------------------------------------------- X* X* Compute the Left and Right Eigenvectors of T X* X* X* Compute the Right eigenvector Matrix: X* X* X NTEST = 9 X RESULT( 9 ) = ULPINV X DO 120 J = 1, N X LWORK( J ) = .TRUE. X 120 CONTINUE X CALL DTREVC( 'R', LWORK, N, T1, LDA, EVECTR, LDU, DUMMA, X $ LDU, N, IN, WORK, IINFO ) X IF( IINFO.NE.0 ) THEN X WRITE( NOUNIT, FMT = 9999 )'DTREVC(R)', IINFO, N, JTYPE, X $ IOLDSD X INFO = ABS( IINFO ) X GO TO 160 X END IF X* X* X* Test 9: | TR - RW | / ( |T| |R| ulp ) X* X CALL DGET22( 'N', 'N', 'N', N, T1, LDA, EVECTR, LDU, WR1, X $ WI1, WORK, DUMMA( 1 ) ) X RESULT( 9 ) = DUMMA( 1 ) X IF( DUMMA( 2 ).GT.THRESH ) THEN X WRITE( NOUNIT, FMT = 9998 )'Right', 'DTREVC', X $ DUMMA( 2 ), N, JTYPE, IOLDSD X END IF X* X* X* X* Compute the Left eigenvector Matrix: X* X* X NTEST = 10 X RESULT( 10 ) = ULPINV X DO 130 J = 1, N X LWORK( J ) = .TRUE. X 130 CONTINUE X CALL DTREVC( 'L', LWORK, N, T1, LDA, DUMMA, LDU, EVECTL, X $ LDU, N, IN, WORK, IINFO ) X IF( IINFO.NE.0 ) THEN X WRITE( NOUNIT, FMT = 9999 )'DTREVC(L)', IINFO, N, JTYPE, X $ IOLDSD X INFO = ABS( IINFO ) X GO TO 160 X END IF X* X* X* Test 10: | LT - WL | / ( |T| |L| ulp ) X* X CALL DGET22( 'Trans', 'N', 'Conj', N, T1, LDA, EVECTL, LDU, X $ WR1, WI1, WORK, DUMMA( 3 ) ) X RESULT( 10 ) = DUMMA( 3 ) X IF( DUMMA( 4 ).GT.THRESH ) THEN X WRITE( NOUNIT, FMT = 9998 )'Left', 'DTREVC', DUMMA( 4 ), X $ N, JTYPE, IOLDSD X END IF X* X*----------------------------------------------------------------------- X* X* X* Call DHSEIN for Right eigenvectors of H, do test 11 X* X NTEST = 11 X RESULT( 11 ) = ULPINV X DO 140 J = 1, N X LWORK( J ) = .TRUE. X 140 CONTINUE X* X CALL DHSEIN( 'R', LWORK, 'N', 'N', N, H, LDA, WR3, WI3, X $ EVECTX, LDU, DUMMA, LDU, N1, IN, WORK, IINFO ) X IF( IINFO.NE.0 ) THEN X WRITE( NOUNIT, FMT = 9999 )'DHSEIN(R)', IINFO, N, JTYPE, X $ IOLDSD X INFO = ABS( IINFO ) X IF( IINFO.LT.0 ) X $ GO TO 160 X ELSE X* X* Test 11: | HX - XW | / ( |H| |X| ulp ) X* X* (from inverse iteration) X* X CALL DGET22( 'N', 'N', 'N', N, H, LDA, EVECTX, LDU, WR3, X $ WI3, WORK, DUMMA( 1 ) ) X IF( DUMMA( 1 ).LT.ULPINV ) X $ RESULT( 11 ) = DUMMA( 1 )*ANINV X IF( DUMMA( 2 ).GT.THRESH ) THEN X WRITE( NOUNIT, FMT = 9998 )'Right', 'DHSEIN', X $ DUMMA( 2 ), N, JTYPE, IOLDSD X END IF X END IF X* X* X* Call DHSEIN for Left eigenvectors of H, do test 12 X* X NTEST = 12 X RESULT( 12 ) = ULPINV X DO 150 J = 1, N X LWORK( J ) = .TRUE. X 150 CONTINUE X* X CALL DHSEIN( 'L', LWORK, 'N', 'N', N, H, LDA, WR3, WI3, X $ DUMMA, LDU, EVECTY, LDU, N1, IN, WORK, IINFO ) X IF( IINFO.NE.0 ) THEN X WRITE( NOUNIT, FMT = 9999 )'DHSEIN(L)', IINFO, N, JTYPE, X $ IOLDSD X INFO = ABS( IINFO ) X IF( IINFO.LT.0 ) X $ GO TO 160 X ELSE X* X* Test 12: | YH - WY | / ( |H| |Y| ulp ) X* X* (from inverse iteration) X* X CALL DGET22( 'C', 'N', 'C', N, H, LDA, EVECTY, LDU, WR3, X $ WI3, WORK, DUMMA( 3 ) ) X IF( DUMMA( 3 ).LT.ULPINV ) X $ RESULT( 12 ) = DUMMA( 3 )*ANINV X IF( DUMMA( 4 ).GT.THRESH ) THEN X WRITE( NOUNIT, FMT = 9998 )'Left', 'DHSEIN', X $ DUMMA( 4 ), N, JTYPE, IOLDSD X END IF X END IF X* X* X*----------------------------------------------------------------------- X* X* End of Loop -- Check for RESULT(j) > THRESH X* X 160 CONTINUE X* X NTESTT = NTESTT + NTEST X CALL DLAFTS( 'DHS', N, N, JTYPE, NTEST, RESULT, IOLDSD, X $ THRESH, NOUNIT, NERRS ) X* X 170 CONTINUE X 180 CONTINUE X* X* Summary X* X CALL DLASUM( 'DHS', NOUNIT, NERRS, NTESTT ) X* X* X*----------------------------------------------------------------------- X* X* X 9999 FORMAT( ' DCHK21: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', X $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) X 9998 FORMAT( ' DCHK21: ', A, ' Eigenvectors from ', A, ' incorrectly ', X $ 'normalized.', / ' Bits of error=', 0P, G10.3, ',', 9X, X $ 'N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, X $ ')' ) X* X RETURN X* X* End of DCHK21 X* X END END_OF_FILE if test 32890 -ne `wc -c <'dchk21.f'`; then echo shar: \"'dchk21.f'\" unpacked with wrong size! fi # end of 'dchk21.f' fi if test -f 'dcopy.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'dcopy.f'\" else echo shar: Extracting \"'dcopy.f'\" \(1490 characters\) sed "s/^X//" >'dcopy.f' <<'END_OF_FILE' X SUBROUTINE DCOPY( N, DX, INCX, DY, INCY ) X* X* copies a vector, x, to a vector, y. X* uses unrolled loops for increments equal to one. X* jack dongarra, linpack, 3/11/78. X* X* .. Scalar Arguments .. X INTEGER INCX, INCY, N X* .. X* .. Array Arguments .. X DOUBLE PRECISION DX( 1 ), DY( 1 ) X* .. X* .. Local Scalars .. X INTEGER I, IX, IY, M, MP1 X* .. X* .. Intrinsic Functions .. X INTRINSIC MOD X* .. X* .. Executable Statements .. X* X IF( N.LE.0 ) X $ RETURN X IF( INCX.EQ.1 .AND. INCY.EQ.1 ) X $ GO TO 20 X* X* code for unequal increments or equal increments X* not equal to 1 X* X IX = 1 X IY = 1 X IF( INCX.LT.0 ) X $ IX = ( -N+1 )*INCX + 1 X IF( INCY.LT.0 ) X $ IY = ( -N+1 )*INCY + 1 X DO 10 I = 1, N X DY( IY ) = DX( IX ) X IX = IX + INCX X IY = IY + INCY X 10 CONTINUE X RETURN X* X* code for both increments equal to 1 X* X* X* clean-up loop X* X 20 M = MOD( N, 7 ) X IF( M.EQ.0 ) X $ GO TO 40 X DO 30 I = 1, M X DY( I ) = DX( I ) X 30 CONTINUE X IF( N.LT.7 ) X $ RETURN X 40 MP1 = M + 1 X DO 50 I = MP1, N, 7 X DY( I ) = DX( I ) X DY( I+1 ) = DX( I+1 ) X DY( I+2 ) = DX( I+2 ) X DY( I+3 ) = DX( I+3 ) X DY( I+4 ) = DX( I+4 ) X DY( I+5 ) = DX( I+5 ) X DY( I+6 ) = DX( I+6 ) X 50 CONTINUE X RETURN X END END_OF_FILE if test 1490 -ne `wc -c <'dcopy.f'`; then echo shar: \"'dcopy.f'\" unpacked with wrong size! fi # end of 'dcopy.f' fi if test -f 'ddot.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'ddot.f'\" else echo shar: Extracting \"'ddot.f'\" \(1663 characters\) sed "s/^X//" >'ddot.f' <<'END_OF_FILE' X DOUBLE PRECISION FUNCTION DDOT( N, DX, INCX, DY, INCY ) X* X* forms the dot product of two vectors. X* uses unrolled loops for increments equal to one. X* jack dongarra, linpack, 3/11/78. X* X* .. Scalar Arguments .. X INTEGER INCX, INCY, N X* .. X* .. Array Arguments .. X DOUBLE PRECISION DX( 1 ), DY( 1 ) X* .. X* .. Local Scalars .. X INTEGER I, IX, IY, M, MP1 X DOUBLE PRECISION DTEMP X* .. X* .. Intrinsic Functions .. X INTRINSIC MOD X* .. X* .. Executable Statements .. X* X DDOT = 0.0D0 X DTEMP = 0.0D0 X IF( N.LE.0 ) X $ RETURN X IF( INCX.EQ.1 .AND. INCY.EQ.1 ) X $ GO TO 20 X* X* code for unequal increments or equal increments X* not equal to 1 X* X IX = 1 X IY = 1 X IF( INCX.LT.0 ) X $ IX = ( -N+1 )*INCX + 1 X IF( INCY.LT.0 ) X $ IY = ( -N+1 )*INCY + 1 X DO 10 I = 1, N X DTEMP = DTEMP + DX( IX )*DY( IY ) X IX = IX + INCX X IY = IY + INCY X 10 CONTINUE X DDOT = DTEMP X RETURN X* X* code for both increments equal to 1 X* X* X* clean-up loop X* X 20 M = MOD( N, 5 ) X IF( M.EQ.0 ) X $ GO TO 40 X DO 30 I = 1, M X DTEMP = DTEMP + DX( I )*DY( I ) X 30 CONTINUE X IF( N.LT.5 ) X $ GO TO 60 X 40 MP1 = M + 1 X DO 50 I = MP1, N, 5 X DTEMP = DTEMP + DX( I )*DY( I ) + DX( I+1 )*DY( I+1 ) + X $ DX( I+2 )*DY( I+2 ) + DX( I+3 )*DY( I+3 ) + X $ DX( I+4 )*DY( I+4 ) X 50 CONTINUE X 60 DDOT = DTEMP X RETURN X END END_OF_FILE if test 1663 -ne `wc -c <'ddot.f'`; then echo shar: \"'ddot.f'\" unpacked with wrong size! fi # end of 'ddot.f' fi if test -f 'depslon.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'depslon.f'\" else echo shar: Extracting \"'depslon.f'\" \(1334 characters\) sed "s/^X//" >'depslon.f' <<'END_OF_FILE' X double precision function epslon (x) X double precision x c c estimate unit roundoff in quantities of size x. c X double precision a,b,c,eps c c this program should function properly on all systems c satisfying the following two assumptions, c 1. the base used in representing floating point c numbers is not a power of three. c 2. the quantity a in statement 10 is represented to c the accuracy used in floating point variables c that are stored in memory. c the statement number 10 and the go to 10 are intended to c force optimizing compilers to generate code satisfying c assumption 2. c under these assumptions, it should be true that, c a is not exactly equal to four-thirds, c b has a zero for its last bit or digit, c c is not exactly equal to one, c eps measures the separation of 1.0 from c the next larger floating point number. c the developers of eispack would appreciate being informed c about any systems where these assumptions do not hold. c c this version dated 4/6/83. c X a = 4.0d0/3.0d0 X 10 b = a - 1.0d0 X c = b + b + b X eps = dabs(c-1.0d0) X if (eps .eq. 0.0d0) go to 10 X epslon = eps*dabs(x) X return X end END_OF_FILE if test 1334 -ne `wc -c <'depslon.f'`; then echo shar: \"'depslon.f'\" unpacked with wrong size! fi # end of 'depslon.f' fi if test -f 'dgehd3.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'dgehd3.f'\" else echo shar: Extracting \"'dgehd3.f'\" \(10674 characters\) sed "s/^X//" >'dgehd3.f' <<'END_OF_FILE' X SUBROUTINE DGEHD3( JOB, N, A, LDA, U, LDU, S, WORK, LDWORK, NWORK, X $ INFO ) X* X* -- LAPACK routine (preliminary version) -- X* Univ. of Tennessee, Oak Ridge National Lab, Argonne National Lab, X* Courant Institute, NAG Ltd., and Rice University X* March 26, 1990 X* X* .. Scalar Arguments .. X CHARACTER JOB X INTEGER INFO, LDA, LDU, LDWORK, N, NWORK X* .. X* X* .. Array Arguments .. X DOUBLE PRECISION A( LDA, * ), S( * ), U( LDU, * ), X $ WORK( LDWORK, * ) X* .. X* X* Purpose X* ======= X* X* This subroutine reduces an N-by-N real general matrix A to X* upper Hessenberg form, and returns the orthogonal transforma- X* tion matrix if desired. X* X* A will be decomposed as U H U', where U is orthogonal, H is X* upper Hessenberg, and U' denotes the transpose of U. X* X* Arguments X* ========= X* X* JOB - CHARACTER*1 X* JOB specifies the computation to be done by DGEHD3. X* JOB = 'H': return the upper Hessenberg matrix H only. X* JOB = 'Q': return both the upper Hessenberg matrix H X* and the orthogonal matrix U. X* Not modified. X* X* N - INTEGER X* N specifies the order of the matrix A. X* N must be at least zero. X* Not modified. X* X* A - DOUBLE PRECISION array, dimension (LDA,N) X* On entry, A specifies the array which contains the matrix X* being reduced. X* On exit, the array A is overwritten by its Hessenberg form. X* X* LDA - INTEGER X* LDA specifies the first dimension of A as X* declared in the calling (sub)program. LDA must be at X* least max(1, N). X* Not modified. X* X* U - DOUBLE PRECISION array, dimension (LDU,N) X* If JOB='Q', then on exit, the N-by-N matrix U will contain X* the orthogonal matrix U used to reduce A to Hessenberg form. X* If JOB = 'H', U is not referenced. X* X* LDU - INTEGER X* LDU specifies the first dimensiion of U as declared in X* the calling (sub)program. LDU must be at least max(1, N). X* If JOB = 'H', U is not referenced. X* Not modified. X* X* S - DOUBLE PRECISION array, dimension (N). X* Workspace. If JOB = 'H', S is not referenced. X* X* WORK - DOUBLE PRECISION array, dimension (LDWORK,NWORK) X* Workspace. X* X* LDWORK - INTEGER X* LDWORK specifies the first dimension of WORK as X* declared in the calling (sub)program. LDWORK must be at X* least max(1, N). X* Not modified. X* X* NWORK - INTEGER X* NWORK specifies the number of columns in WORK. X* NWORK must be at least 4, and should be at least 3*NB+1, X* where NB is the blocksize as returned by the routine ENVIR. X* Not modified. X* X* INFO - INTEGER X* On exit, INFO is set to X* 0 normal return. X* -k if input argument number k is illegal. X* X* Internal parameters which may be modified by the user X* ===================================================== X* X* NB - INTEGER X* NB specifies the block size. It is normally gotten X* by a call to the subroutine ENVIR. X* Not modified. X* X* .. Parameters .. X DOUBLE PRECISION ZERO, ONE X PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) X* .. X* X* .. Local Scalars .. X INTEGER I, IFST, IJOB, ILST, IRES, J, JIFST, JK, JS, X $ JS1, KB, LS, LS3, M, NB, NBL, NN, NNB X DOUBLE PRECISION DELTA, TAU X* .. X* X* .. External Functions .. X LOGICAL LSAME X DOUBLE PRECISION DDOT X EXTERNAL LSAME, DDOT X* .. X* X* .. External Subroutines .. X EXTERNAL ENVIR, DAXPY, DGEMM, DGEMV, DLARFG, DLAZRO, X $ DORGC3, DSCAL, XERBLA X* .. X* X* .. Intrinsic Functions .. X INTRINSIC MAX, MIN, MOD X* .. X* X* .. Executable Statements .. X* X* See "Block Reduction of Matrices to condensed Forms for X* Eigenvalue Computation" by J. Dongarra, S. Hammarling and X* D. Sorensen, LAPACK Working Note #2, and "On a Block X* Implementation of the Hessenberg Multishift QR Iteration" X* by Z. Bai and J. Demmel, LAPACK Working Note #8 for a X* detailed description of the algorithm. X* X* Decode and Test the input parameters X* X IF( LSAME( JOB, 'H' ) ) THEN X IJOB = 1 X ELSE IF( LSAME( JOB, 'Q' ) ) THEN X IJOB = 2 X ELSE X IJOB = -1 X END IF X* X INFO = 0 X IF( IJOB.EQ.-1 ) THEN X INFO = -1 X ELSE IF( N.LT.0 ) THEN X INFO = -2 X ELSE IF( LDA.LT.MAX( 1, N ) ) THEN X INFO = -4 X ELSE IF( LDWORK.LT.MAX( 1, N ) ) THEN X INFO = -8 X END IF X IF( IJOB.EQ.2 ) THEN X IF( LDU.LT.MAX( 1, N ) ) X $ INFO = -6 X END IF X IF( INFO.NE.0 ) THEN X CALL XERBLA( 'DGEHD3', -INFO ) X RETURN X END IF X* X* Initialize U, if desired. X* X IF( IJOB.EQ.2 ) THEN X CALL DLAZRO( N, N, ZERO, ZERO, U, LDU ) X END IF X* X* Quick return if possible. X* X IF( N.LE.2 ) X $ GO TO 130 X* X* Get block sizes. X* X CALL ENVIR( 'BLOCK', NB ) X NBL = MIN( NB, N-2, ( NWORK-1 ) / 3 ) X* X* Determine the number of blocks: NNB X* X IRES = MOD( N-2, NBL ) X NNB = ( N-2 ) / NBL X NN = NNB X IF( IRES.NE.0 ) X $ NNB = NNB + 1 X* X* Outer loop in block updating X* X DO 100 KB = 1, NNB X* X* Determine the first column index IFST and the X* last column index ILST of the KBth block. X* X IFST = ( KB-1 )*NBL + 1 X ILST = IFST + NBL - 1 X IF( KB.EQ.NN+1 ) X $ ILST = N - 2 X LS = ILST - IFST + 1 X LS3 = 3*LS + 1 X* X* Initialize working array WORK. Note that in the working X* array WORK: WORK(*,1:LS) store V; WORK(*,LS+1:2*LS) X* store W and WORK(*,1+2*LS:3*LS) store U as described in X* the above papers. On each block, matrix A will be X* updated as X* A = A - U*V' - W*U' X* X DO 40 J = 1, LS X DO 10 I = IFST, N X WORK( I, J ) = ZERO X 10 CONTINUE X DO 20 I = 1, N X WORK( I, J+LS ) = ZERO X 20 CONTINUE X DO 30 I = IFST, N X WORK( I, J+2*LS ) = ZERO X 30 CONTINUE X 40 CONTINUE X* X* Inner loop in each block X* X DO 90 JK = IFST, ILST X* X* form the JKth column. X* X JIFST = JK - IFST + 1 + 2*LS X DO 50 I = JK + 1, N X WORK( I, JIFST ) = A( I, JK ) X 50 CONTINUE X* X DO 70 J = IFST, JK - 1 X DO 60 I = JK + 1, N X WORK( I, JIFST ) = WORK( I, JIFST ) - X $ WORK( JK, J-IFST+1 )* X $ WORK( I, J-IFST+1+2*LS ) - X $ WORK( JK, J-IFST+1+2*LS )* X $ WORK( I, J-IFST+1+LS ) X 60 CONTINUE X 70 CONTINUE X* X* Compute Householder transformation for column JK: X* X CALL DLARFG( N-JK, WORK( JK+1, JIFST ), WORK( JK+2, JIFST ), X $ 1, TAU ) X WORK( JK+1, JIFST ) = ONE X* X IF( IJOB.EQ.2 ) THEN X DO 80 I = JK + 1, N X U( I, JK ) = WORK( I, JIFST ) X 80 CONTINUE X S( JK ) = TAU X END IF X* X* Aggregate the transformation vectors in inner loop. X* X* A'*uj - V*U'*uj - U*W'*uj --> vj X* X JS = JK - IFST X JS1 = JS + 1 X CALL DGEMV( 'T', N-JK, JS, ONE, WORK( JK+1, 1+LS ), LDWORK, X $ WORK( JK+1, JIFST ), 1, ZERO, WORK( 1, LS3 ), X $ 1 ) X* X CALL DGEMV( 'N', N-IFST, JS, ONE, WORK( IFST+1, 1+2*LS ), X $ LDWORK, WORK( 1, LS3 ), 1, ZERO, X $ WORK( IFST+1, JS1 ), 1 ) X* X CALL DGEMV( 'T', N-JK, JS, ONE, WORK( JK+1, 1+2*LS ), X $ LDWORK, WORK( JK+1, JIFST ), 1, ZERO, X $ WORK( 1, LS3 ), 1 ) X* X CALL DGEMV( 'N', N-IFST+1, JS, ONE, WORK( IFST, 1 ), LDWORK, X $ WORK( 1, LS3 ), 1, ONE, WORK( IFST, JS1 ), 1 ) X* X CALL DGEMV( 'T', N-JK, N-IFST+1, ONE, A( JK+1, IFST ), LDA, X $ WORK( JK+1, JIFST ), 1, -ONE, WORK( IFST, JS1 ), X $ 1 ) X* X* A*uj - U*V'*uj - W*U'*uj --> wj X* X CALL DGEMV( 'N', N, JS, ONE, WORK( 1, LS+1 ), LDWORK, X $ WORK( 1, LS3 ), 1, ZERO, WORK( 1, JS1+LS ), 1 ) X* X CALL DGEMV( 'T', N-JK, JS, ONE, WORK( JK+1, 1 ), LDWORK, X $ WORK( JK+1, JIFST ), 1, ZERO, WORK( 1, LS3 ), X $ 1 ) X* X CALL DGEMV( 'N', N-IFST, JS, ONE, WORK( IFST+1, 1+2*LS ), X $ LDWORK, WORK( 1, LS3 ), 1, ONE, X $ WORK( IFST+1, JS1+LS ), 1 ) X* X CALL DGEMV( 'N', N, N-JK, ONE, A( 1, JK+1 ), LDA, X $ WORK( JK+1, JIFST ), 1, -ONE, WORK( 1, JS1+LS ), X $ 1 ) X* X DELTA = DDOT( N-JK, WORK( JK+1, JS1+LS ), 1, X $ WORK( JK+1, JIFST ), 1 ) X CALL DAXPY( N-JK, -TAU*DELTA, WORK( JK+1, JIFST ), 1, X $ WORK( JK+1, JS1 ), 1 ) X* X CALL DSCAL( N-IFST+1, TAU, WORK( IFST, JS1 ), 1 ) X CALL DSCAL( N, TAU, WORK( 1, JS1+LS ), 1 ) X* X 90 CONTINUE X* X* The end of inner JK loop X* X* Row block updating: X* A = A(IFST+1:N,IFST:N) - U*V' X* X CALL DGEMM( 'N', 'T', N-IFST, N-IFST+1, LS, -ONE, X $ WORK( IFST+1, 1+2*LS ), LDWORK, WORK( IFST, 1 ), X $ LDWORK, ONE, A( IFST+1, IFST ), LDA ) X* X* Column block updating: X* A = A(1:N,IFST+1:N) - W*U' X* X CALL DGEMM( 'N', 'T', N, N-IFST, LS, -ONE, WORK( 1, 1+LS ), X $ LDWORK, WORK( IFST+1, 1+2*LS ), LDWORK, ONE, X $ A( 1, IFST+1 ), LDA ) X* X 100 CONTINUE X* X* Clean up X* X DO 120 J = 1, N - 2 X DO 110 I = J + 2, N X A( I, J ) = ZERO X 110 CONTINUE X 120 CONTINUE X* X* Form orthogonal transformation if desired. X* X 130 CONTINUE X IF( IJOB.EQ.2 ) THEN X M = MAX( N-2, 0 ) X CALL DORGC3( N, M, U, LDU, S, WORK, INFO ) X END IF X* X RETURN X* X* End of DGEHD3 X* X END END_OF_FILE if test 10674 -ne `wc -c <'dgehd3.f'`; then echo shar: \"'dgehd3.f'\" unpacked with wrong size! fi # end of 'dgehd3.f' fi if test -f 'dgemm.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'dgemm.f'\" else echo shar: Extracting \"'dgemm.f'\" \(9851 characters\) sed "s/^X//" >'dgemm.f' <<'END_OF_FILE' X SUBROUTINE DGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, X $ BETA, C, LDC ) X* .. Scalar Arguments .. X CHARACTER*1 TRANSA, TRANSB X INTEGER M, N, K, LDA, LDB, LDC X DOUBLE PRECISION ALPHA, BETA X* .. Array Arguments .. X DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ) X* .. X* X* Purpose X* ======= X* X* DGEMM performs one of the matrix-matrix operations X* X* C := alpha*op( A )*op( B ) + beta*C, X* X* where op( X ) is one of X* X* op( X ) = X or op( X ) = X', X* X* alpha and beta are scalars, and A, B and C are matrices, with op( A ) X* an m by k matrix, op( B ) a k by n matrix and C an m by n matrix. X* X* Parameters X* ========== X* X* TRANSA - CHARACTER*1. X* On entry, TRANSA specifies the form of op( A ) to be used in X* the matrix multiplication as follows: X* X* TRANSA = 'N' or 'n', op( A ) = A. X* X* TRANSA = 'T' or 't', op( A ) = A'. X* X* TRANSA = 'C' or 'c', op( A ) = A'. X* X* Unchanged on exit. X* X* TRANSB - CHARACTER*1. X* On entry, TRANSB specifies the form of op( B ) to be used in X* the matrix multiplication as follows: X* X* TRANSB = 'N' or 'n', op( B ) = B. X* X* TRANSB = 'T' or 't', op( B ) = B'. X* X* TRANSB = 'C' or 'c', op( B ) = B'. X* X* Unchanged on exit. X* X* M - INTEGER. X* On entry, M specifies the number of rows of the matrix X* op( A ) and of the matrix C. M must be at least zero. X* Unchanged on exit. X* X* N - INTEGER. X* On entry, N specifies the number of columns of the matrix X* op( B ) and the number of columns of the matrix C. N must be X* at least zero. X* Unchanged on exit. X* X* K - INTEGER. X* On entry, K specifies the number of columns of the matrix X* op( A ) and the number of rows of the matrix op( B ). K must X* be at least zero. X* Unchanged on exit. X* X* ALPHA - DOUBLE PRECISION. X* On entry, ALPHA specifies the scalar alpha. X* Unchanged on exit. X* X* A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is X* k when TRANSA = 'N' or 'n', and is m otherwise. X* Before entry with TRANSA = 'N' or 'n', the leading m by k X* part of the array A must contain the matrix A, otherwise X* the leading k by m part of the array A must contain the X* matrix A. X* Unchanged on exit. X* X* LDA - INTEGER. X* On entry, LDA specifies the first dimension of A as declared X* in the calling (sub) program. When TRANSA = 'N' or 'n' then X* LDA must be at least max( 1, m ), otherwise LDA must be at X* least max( 1, k ). X* Unchanged on exit. X* X* B - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is X* n when TRANSB = 'N' or 'n', and is k otherwise. X* Before entry with TRANSB = 'N' or 'n', the leading k by n X* part of the array B must contain the matrix B, otherwise X* the leading n by k part of the array B must contain the X* matrix B. X* Unchanged on exit. X* X* LDB - INTEGER. X* On entry, LDB specifies the first dimension of B as declared X* in the calling (sub) program. When TRANSB = 'N' or 'n' then X* LDB must be at least max( 1, k ), otherwise LDB must be at X* least max( 1, n ). X* Unchanged on exit. X* X* BETA - DOUBLE PRECISION. X* On entry, BETA specifies the scalar beta. When BETA is X* supplied as zero then C need not be set on input. X* Unchanged on exit. X* X* C - DOUBLE PRECISION array of DIMENSION ( LDC, n ). X* Before entry, the leading m by n part of the array C must X* contain the matrix C, except when beta is zero, in which X* case C need not be set on entry. X* On exit, the array C is overwritten by the m by n matrix X* ( alpha*op( A )*op( B ) + beta*C ). X* X* LDC - INTEGER. X* On entry, LDC specifies the first dimension of C as declared X* in the calling (sub) program. LDC must be at least X* max( 1, m ). X* Unchanged on exit. X* X* X* Level 3 Blas routine. X* X* -- Written on 8-February-1989. X* Jack Dongarra, Argonne National Laboratory. X* Iain Duff, AERE Harwell. X* Jeremy Du Croz, Numerical Algorithms Group Ltd. X* Sven Hammarling, Numerical Algorithms Group Ltd. X* X* X* .. External Functions .. X LOGICAL LSAME X EXTERNAL LSAME X* .. External Subroutines .. X EXTERNAL XERBLA X* .. Intrinsic Functions .. X INTRINSIC MAX X* .. Local Scalars .. X LOGICAL NOTA, NOTB X INTEGER I, INFO, J, L, NCOLA, NROWA, NROWB X DOUBLE PRECISION TEMP X* .. Parameters .. X DOUBLE PRECISION ONE , ZERO X PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) X* .. X* .. Executable Statements .. X* X* Set NOTA and NOTB as true if A and B respectively are not X* transposed and set NROWA, NCOLA and NROWB as the number of rows X* and columns of A and the number of rows of B respectively. X* X NOTA = LSAME( TRANSA, 'N' ) X NOTB = LSAME( TRANSB, 'N' ) X IF( NOTA )THEN X NROWA = M X NCOLA = K X ELSE X NROWA = K X NCOLA = M X END IF X IF( NOTB )THEN X NROWB = K X ELSE X NROWB = N X END IF X* X* Test the input parameters. X* X INFO = 0 X IF( ( .NOT.NOTA ).AND. X $ ( .NOT.LSAME( TRANSA, 'C' ) ).AND. X $ ( .NOT.LSAME( TRANSA, 'T' ) ) )THEN X INFO = 1 X ELSE IF( ( .NOT.NOTB ).AND. X $ ( .NOT.LSAME( TRANSB, 'C' ) ).AND. X $ ( .NOT.LSAME( TRANSB, 'T' ) ) )THEN X INFO = 2 X ELSE IF( M .LT.0 )THEN X INFO = 3 X ELSE IF( N .LT.0 )THEN X INFO = 4 X ELSE IF( K .LT.0 )THEN X INFO = 5 X ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN X INFO = 8 X ELSE IF( LDB.LT.MAX( 1, NROWB ) )THEN X INFO = 10 X ELSE IF( LDC.LT.MAX( 1, M ) )THEN X INFO = 13 X END IF X IF( INFO.NE.0 )THEN X CALL XERBLA( 'DGEMM ', INFO ) X RETURN X END IF X* X* Quick return if possible. X* X IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. X $ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) ) X $ RETURN X* X* And if alpha.eq.zero. X* X IF( ALPHA.EQ.ZERO )THEN X IF( BETA.EQ.ZERO )THEN X DO 20, J = 1, N X DO 10, I = 1, M X C( I, J ) = ZERO X 10 CONTINUE X 20 CONTINUE X ELSE X DO 40, J = 1, N X DO 30, I = 1, M X C( I, J ) = BETA*C( I, J ) X 30 CONTINUE X 40 CONTINUE X END IF X RETURN X END IF X* X* Start the operations. X* X IF( NOTB )THEN X IF( NOTA )THEN X* X* Form C := alpha*A*B + beta*C. X* X DO 90, J = 1, N X IF( BETA.EQ.ZERO )THEN X DO 50, I = 1, M X C( I, J ) = ZERO X 50 CONTINUE X ELSE IF( BETA.NE.ONE )THEN X DO 60, I = 1, M X C( I, J ) = BETA*C( I, J ) X 60 CONTINUE X END IF X DO 80, L = 1, K X IF( B( L, J ).NE.ZERO )THEN X TEMP = ALPHA*B( L, J ) X DO 70, I = 1, M X C( I, J ) = C( I, J ) + TEMP*A( I, L ) X 70 CONTINUE X END IF X 80 CONTINUE X 90 CONTINUE X ELSE X* X* Form C := alpha*A'*B + beta*C X* X DO 120, J = 1, N X DO 110, I = 1, M X TEMP = ZERO X DO 100, L = 1, K X TEMP = TEMP + A( L, I )*B( L, J ) X 100 CONTINUE X IF( BETA.EQ.ZERO )THEN X C( I, J ) = ALPHA*TEMP X ELSE X C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) X END IF X 110 CONTINUE X 120 CONTINUE X END IF X ELSE X IF( NOTA )THEN X* X* Form C := alpha*A*B' + beta*C X* X DO 170, J = 1, N X IF( BETA.EQ.ZERO )THEN X DO 130, I = 1, M X C( I, J ) = ZERO X 130 CONTINUE X ELSE IF( BETA.NE.ONE )THEN X DO 140, I = 1, M X C( I, J ) = BETA*C( I, J ) X 140 CONTINUE X END IF X DO 160, L = 1, K X IF( B( J, L ).NE.ZERO )THEN X TEMP = ALPHA*B( J, L ) X DO 150, I = 1, M X C( I, J ) = C( I, J ) + TEMP*A( I, L ) X 150 CONTINUE X END IF X 160 CONTINUE X 170 CONTINUE X ELSE X* X* Form C := alpha*A'*B' + beta*C X* X DO 200, J = 1, N X DO 190, I = 1, M X TEMP = ZERO X DO 180, L = 1, K X TEMP = TEMP + A( L, I )*B( J, L ) X 180 CONTINUE X IF( BETA.EQ.ZERO )THEN X C( I, J ) = ALPHA*TEMP X ELSE X C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) X END IF X 190 CONTINUE X 200 CONTINUE X END IF X END IF X* X RETURN X* X* End of DGEMM . X* X END END_OF_FILE if test 9851 -ne `wc -c <'dgemm.f'`; then echo shar: \"'dgemm.f'\" unpacked with wrong size! fi # end of 'dgemm.f' fi if test -f 'dgemv.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'dgemv.f'\" else echo shar: Extracting \"'dgemv.f'\" \(7481 characters\) sed "s/^X//" >'dgemv.f' <<'END_OF_FILE' X SUBROUTINE DGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX, X $ BETA, Y, INCY ) X* .. Scalar Arguments .. X DOUBLE PRECISION ALPHA, BETA X INTEGER INCX, INCY, LDA, M, N X CHARACTER*1 TRANS X* .. Array Arguments .. X DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) X* .. X* X* Purpose X* ======= X* X* DGEMV performs one of the matrix-vector operations X* X* y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, X* X* where alpha and beta are scalars, x and y are vectors and A is an X* m by n matrix. X* X* Parameters X* ========== X* X* TRANS - CHARACTER*1. X* On entry, TRANS specifies the operation to be performed as X* follows: X* X* TRANS = 'N' or 'n' y := alpha*A*x + beta*y. X* X* TRANS = 'T' or 't' y := alpha*A'*x + beta*y. X* X* TRANS = 'C' or 'c' y := alpha*A'*x + beta*y. X* X* Unchanged on exit. X* X* M - INTEGER. X* On entry, M specifies the number of rows of the matrix A. X* M must be at least zero. X* Unchanged on exit. X* X* N - INTEGER. X* On entry, N specifies the number of columns of the matrix A. X* N must be at least zero. X* Unchanged on exit. X* X* ALPHA - DOUBLE PRECISION. X* On entry, ALPHA specifies the scalar alpha. X* Unchanged on exit. X* X* A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). X* Before entry, the leading m by n part of the array A must X* contain the matrix of coefficients. X* Unchanged on exit. X* X* LDA - INTEGER. X* On entry, LDA specifies the first dimension of A as declared X* in the calling (sub) program. LDA must be at least X* max( 1, m ). X* Unchanged on exit. X* X* X - DOUBLE PRECISION array of DIMENSION at least X* ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' X* and at least X* ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. X* Before entry, the incremented array X must contain the X* vector x. X* Unchanged on exit. X* X* INCX - INTEGER. X* On entry, INCX specifies the increment for the elements of X* X. INCX must not be zero. X* Unchanged on exit. X* X* BETA - DOUBLE PRECISION. X* On entry, BETA specifies the scalar beta. When BETA is X* supplied as zero then Y need not be set on input. X* Unchanged on exit. X* X* Y - DOUBLE PRECISION array of DIMENSION at least X* ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' X* and at least X* ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. X* Before entry with BETA non-zero, the incremented array Y X* must contain the vector y. On exit, Y is overwritten by the X* updated vector y. X* X* INCY - INTEGER. X* On entry, INCY specifies the increment for the elements of X* Y. INCY must not be zero. X* Unchanged on exit. X* X* X* Level 2 Blas routine. X* X* -- Written on 22-October-1986. X* Jack Dongarra, Argonne National Lab. X* Jeremy Du Croz, Nag Central Office. X* Sven Hammarling, Nag Central Office. X* Richard Hanson, Sandia National Labs. X* X* X* .. Parameters .. X DOUBLE PRECISION ONE , ZERO X PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) X* .. Local Scalars .. X DOUBLE PRECISION TEMP X INTEGER I, INFO, IX, IY, J, JX, JY, KX, KY, LENX, LENY X* .. External Functions .. X LOGICAL LSAME X EXTERNAL LSAME X* .. External Subroutines .. X EXTERNAL XERBLA X* .. Intrinsic Functions .. X INTRINSIC MAX X* .. X* .. Executable Statements .. X* X* Test the input parameters. X* X INFO = 0 X IF ( .NOT.LSAME( TRANS, 'N' ).AND. X $ .NOT.LSAME( TRANS, 'T' ).AND. X $ .NOT.LSAME( TRANS, 'C' ) )THEN X INFO = 1 X ELSE IF( M.LT.0 )THEN X INFO = 2 X ELSE IF( N.LT.0 )THEN X INFO = 3 X ELSE IF( LDA.LT.MAX( 1, M ) )THEN X INFO = 6 X ELSE IF( INCX.EQ.0 )THEN X INFO = 8 X ELSE IF( INCY.EQ.0 )THEN X INFO = 11 X END IF X IF( INFO.NE.0 )THEN X CALL XERBLA( 'DGEMV ', INFO ) X RETURN X END IF X* X* Quick return if possible. X* X IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. X $ ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) ) X $ RETURN X* X* Set LENX and LENY, the lengths of the vectors x and y, and set X* up the start points in X and Y. X* X IF( LSAME( TRANS, 'N' ) )THEN X LENX = N X LENY = M X ELSE X LENX = M X LENY = N X END IF X IF( INCX.GT.0 )THEN X KX = 1 X ELSE X KX = 1 - ( LENX - 1 )*INCX X END IF X IF( INCY.GT.0 )THEN X KY = 1 X ELSE X KY = 1 - ( LENY - 1 )*INCY X END IF X* X* Start the operations. In this version the elements of A are X* accessed sequentially with one pass through A. X* X* First form y := beta*y. X* X IF( BETA.NE.ONE )THEN X IF( INCY.EQ.1 )THEN X IF( BETA.EQ.ZERO )THEN X DO 10, I = 1, LENY X Y( I ) = ZERO X 10 CONTINUE X ELSE X DO 20, I = 1, LENY X Y( I ) = BETA*Y( I ) X 20 CONTINUE X END IF X ELSE X IY = KY X IF( BETA.EQ.ZERO )THEN X DO 30, I = 1, LENY X Y( IY ) = ZERO X IY = IY + INCY X 30 CONTINUE X ELSE X DO 40, I = 1, LENY X Y( IY ) = BETA*Y( IY ) X IY = IY + INCY X 40 CONTINUE X END IF X END IF X END IF X IF( ALPHA.EQ.ZERO ) X $ RETURN X IF( LSAME( TRANS, 'N' ) )THEN X* X* Form y := alpha*A*x + y. X* X JX = KX X IF( INCY.EQ.1 )THEN X DO 60, J = 1, N X IF( X( JX ).NE.ZERO )THEN X TEMP = ALPHA*X( JX ) X DO 50, I = 1, M X Y( I ) = Y( I ) + TEMP*A( I, J ) X 50 CONTINUE X END IF X JX = JX + INCX X 60 CONTINUE X ELSE X DO 80, J = 1, N X IF( X( JX ).NE.ZERO )THEN X TEMP = ALPHA*X( JX ) X IY = KY X DO 70, I = 1, M X Y( IY ) = Y( IY ) + TEMP*A( I, J ) X IY = IY + INCY X 70 CONTINUE X END IF X JX = JX + INCX X 80 CONTINUE X END IF X ELSE X* X* Form y := alpha*A'*x + y. X* X JY = KY X IF( INCX.EQ.1 )THEN X DO 100, J = 1, N X TEMP = ZERO X DO 90, I = 1, M X TEMP = TEMP + A( I, J )*X( I ) X 90 CONTINUE X Y( JY ) = Y( JY ) + ALPHA*TEMP X JY = JY + INCY X 100 CONTINUE X ELSE X DO 120, J = 1, N X TEMP = ZERO X IX = KX X DO 110, I = 1, M X TEMP = TEMP + A( I, J )*X( IX ) X IX = IX + INCX X 110 CONTINUE X Y( JY ) = Y( JY ) + ALPHA*TEMP X JY = JY + INCY X 120 CONTINUE X END IF X END IF X* X RETURN X* X* End of DGEMV . X* X END END_OF_FILE if test 7481 -ne `wc -c <'dgemv.f'`; then echo shar: \"'dgemv.f'\" unpacked with wrong size! fi # end of 'dgemv.f' fi if test -f 'dger.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'dger.f'\" else echo shar: Extracting \"'dger.f'\" \(4366 characters\) sed "s/^X//" >'dger.f' <<'END_OF_FILE' X SUBROUTINE DGER ( M, N, ALPHA, X, INCX, Y, INCY, A, LDA ) X* .. Scalar Arguments .. X DOUBLE PRECISION ALPHA X INTEGER INCX, INCY, LDA, M, N X* .. Array Arguments .. X DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) X* .. X* X* Purpose X* ======= X* X* DGER performs the rank 1 operation X* X* A := alpha*x*y' + A, X* X* where alpha is a scalar, x is an m element vector, y is an n element X* vector and A is an m by n matrix. X* X* Parameters X* ========== X* X* M - INTEGER. X* On entry, M specifies the number of rows of the matrix A. X* M must be at least zero. X* Unchanged on exit. X* X* N - INTEGER. X* On entry, N specifies the number of columns of the matrix A. X* N must be at least zero. X* Unchanged on exit. X* X* ALPHA - DOUBLE PRECISION. X* On entry, ALPHA specifies the scalar alpha. X* Unchanged on exit. X* X* X - DOUBLE PRECISION array of dimension at least X* ( 1 + ( m - 1 )*abs( INCX ) ). X* Before entry, the incremented array X must contain the m X* element vector x. X* Unchanged on exit. X* X* INCX - INTEGER. X* On entry, INCX specifies the increment for the elements of X* X. INCX must not be zero. X* Unchanged on exit. X* X* Y - DOUBLE PRECISION array of dimension at least X* ( 1 + ( n - 1 )*abs( INCY ) ). X* Before entry, the incremented array Y must contain the n X* element vector y. X* Unchanged on exit. X* X* INCY - INTEGER. X* On entry, INCY specifies the increment for the elements of X* Y. INCY must not be zero. X* Unchanged on exit. X* X* A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). X* Before entry, the leading m by n part of the array A must X* contain the matrix of coefficients. On exit, A is X* overwritten by the updated matrix. X* X* LDA - INTEGER. X* On entry, LDA specifies the first dimension of A as declared X* in the calling (sub) program. LDA must be at least X* max( 1, m ). X* Unchanged on exit. X* X* X* Level 2 Blas routine. X* X* -- Written on 22-October-1986. X* Jack Dongarra, Argonne National Lab. X* Jeremy Du Croz, Nag Central Office. X* Sven Hammarling, Nag Central Office. X* Richard Hanson, Sandia National Labs. X* X* X* .. Parameters .. X DOUBLE PRECISION ZERO X PARAMETER ( ZERO = 0.0D+0 ) X* .. Local Scalars .. X DOUBLE PRECISION TEMP X INTEGER I, INFO, IX, J, JY, KX X* .. External Subroutines .. X EXTERNAL XERBLA X* .. Intrinsic Functions .. X INTRINSIC MAX X* .. X* .. Executable Statements .. X* X* Test the input parameters. X* X INFO = 0 X IF ( M.LT.0 )THEN X INFO = 1 X ELSE IF( N.LT.0 )THEN X INFO = 2 X ELSE IF( INCX.EQ.0 )THEN X INFO = 5 X ELSE IF( INCY.EQ.0 )THEN X INFO = 7 X ELSE IF( LDA.LT.MAX( 1, M ) )THEN X INFO = 9 X END IF X IF( INFO.NE.0 )THEN X CALL XERBLA( 'DGER ', INFO ) X RETURN X END IF X* X* Quick return if possible. X* X IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) ) X $ RETURN X* X* Start the operations. In this version the elements of A are X* accessed sequentially with one pass through A. X* X IF( INCY.GT.0 )THEN X JY = 1 X ELSE X JY = 1 - ( N - 1 )*INCY X END IF X IF( INCX.EQ.1 )THEN X DO 20, J = 1, N X IF( Y( JY ).NE.ZERO )THEN X TEMP = ALPHA*Y( JY ) X DO 10, I = 1, M X A( I, J ) = A( I, J ) + X( I )*TEMP X 10 CONTINUE X END IF X JY = JY + INCY X 20 CONTINUE X ELSE X IF( INCX.GT.0 )THEN X KX = 1 X ELSE X KX = 1 - ( M - 1 )*INCX X END IF X DO 40, J = 1, N X IF( Y( JY ).NE.ZERO )THEN X TEMP = ALPHA*Y( JY ) X IX = KX X DO 30, I = 1, M X A( I, J ) = A( I, J ) + X( IX )*TEMP X IX = IX + INCX X 30 CONTINUE X END IF X JY = JY + INCY X 40 CONTINUE X END IF X* X RETURN X* X* End of DGER . X* X END END_OF_FILE if test 4366 -ne `wc -c <'dger.f'`; then echo shar: \"'dger.f'\" unpacked with wrong size! fi # end of 'dger.f' fi if test -f 'dget21.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'dget21.f'\" else echo shar: Extracting \"'dget21.f'\" \(8295 characters\) sed "s/^X//" >'dget21.f' <<'END_OF_FILE' X SUBROUTINE DGET21( ITYPE, N, A, LDA, B, LDB, U, LDU, TAU, WORK, X $ RESULT ) X* X* -- LAPACK test routine (preliminary version) -- X* Univ. of Tennessee, Oak Ridge National Lab, Argonne National Lab, X* Courant Institute, NAG Ltd., and Rice University X* March 26, 1990 X* X* .. Scalar Arguments .. X* X INTEGER ITYPE, LDA, LDB, LDU, N X* .. X* X* .. Array Arguments .. X* X DOUBLE PRECISION A( LDA, * ), B( LDB, * ), RESULT( * ), X $ TAU( * ), U( LDU, * ), WORK( * ) X* .. X* X*----------------------------------------------------------------------- X* X* Purpose X* ======= X* X* DGET21 generally checks a decomposition of the form X* X* A = U B U' X* X* where ' means transpose and U is orthogonal. If ITYPE=1, X* then U is represented as a dense matrix, otherwise the X* U is expressed as a product of Householder transformations, X* whose vectors are stored in the array "U" and whose scaling X* constants are in "TAU". X* X* Specifically, if ITYPE=1 or 3, then: X* X* RESULT(1) = | A - U B U' | / ( |A| n ulp ) X* X* If ITYPE=2, then: X* X* RESULT(1) = | A - B | / ( |A| n ulp ) X* X* If ITYPE=4, then: X* X* RESULT(1) = | I - BU' | / ( n ulp ) X* X* If ITYPE=1, then a second check is performed: X* X* RESULT(2) = | I - UU' | / ( n ulp ) X* X* otherwise, RESULT(2) is not modified. X* X* X* For ITYPE > 1, the transformation U is expressed as a product X* U = H(1)...H(n-2), where H(j) = I - tau(j) v(j) v(j)' X* and each vector v(j) has its first j elements 0, the (j+1)st X* assumed to be 1, and the remaining n-1-j elements stored in X* U(n-1-j:n,j). X* X* Arguments X* ========== X* X* ITYPE - INTEGER X* Specifies the type of tests to be performed. X* 1: U expressed as a dense orthogonal matrix: X* RESULT(1) = | A - U B U' | / ( |A| n ulp ) *and* X* RESULT(2) = | I - UU' | / ( n ulp ) X* X* 2: RESULT(1) = | A - B | / ( |A| n ulp ) X* X* 3: U expressed as a product of Housholder transformations: X* RESULT(1) = | A - U B U' | / ( |A| n ulp ) X* X* 4: U expressed as a product of Housholder transformations: X* RESULT(1) = | I - BU' | / ( n ulp ) X* X* N - INTEGER X* The size of the matrix. If it is zero, DGET21 does nothing. X* It must be at least zero. X* Not modified. X* X* A - DOUBLE PRECISION array of dimension ( LDA , N ) X* The original (unfactored) matrix. X* Not referenced if ITYPE=4. X* Not modified. X* X* LDA - INTEGER X* The leading dimension of A. It must be at least 1 X* and at least N. X* Not modified. X* X* B - DOUBLE PRECISION array of dimension ( LDB , N ) X* The factored matrix. X* Not modified. X* X* LDB - INTEGER X* The leading dimension of B. It must be at least 1 X* and at least N. X* Not modified. X* X* U - DOUBLE PRECISION array of dimension ( LDU, N ). X* The orthogonal matrix in the decomposition. If ITYPE=1, X* then it is just the matrix, otherwise the lower triangle X* contains the Householder vectors used to describe U. X* Not referenced if ITYPE=2 X* Not modified. X* X* LDU - INTEGER X* The leading dimension of U. LDU must be at least N and X* at least 1. X* Not modified. X* X* TAU - DOUBLE PRECISION array of dimension ( N ) X* If ITYPE > 2, then TAU(j) is the scalar factor of X* v(j) v(j)' in the Householder transformation H(j) of X* the product U = H(1)...H(n-2) X* If ITYPE <= 2, then TAU is not referenced. X* Not modified. X* X* WORK - DOUBLE PRECISION array of dimension ( 2*N**2 ) X* Workspace. X* Modified. X* X* RESULT - DOUBLE PRECISION array of dimension ( 2 ) X* The values computed by the two tests described above. The X* values are currently limited to 1/ulp, to avoid overflow. X* Errors are flagged by RESULT(1)=10/ulp. X* RESULT(1) is always modified. RESULT(2) is modified only X* if ITYPE=1. X* Modified. X* X*----------------------------------------------------------------------- X* X* X* .. Parameters .. X* X DOUBLE PRECISION ZERO, ONE, TEN X PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TEN = 10.0D0 ) X* .. X* X* .. Local Scalars .. X* X INTEGER IINFO, JCOL, JDIAG, JROW X DOUBLE PRECISION ANORM, ULP, UNFL, WNORM X* .. X* X* .. External Functions .. X* X DOUBLE PRECISION DLAMCH, DLANGE X EXTERNAL DLAMCH, DLANGE X* .. X* X* .. External Subroutines .. X* X EXTERNAL DGEMM, DLACPY, DORMC2 X* .. X* X* .. Intrinsic Functions .. X* X INTRINSIC DBLE, MAX, MIN X* .. X* X* X*----------------------------------------------------------------------- X* .. Executable Statements .. X* X RESULT( 1 ) = ZERO X IF( ITYPE.EQ.1 ) X $ RESULT( 2 ) = ZERO X IF( N.LE.0 ) X $ RETURN X* X* Constants X* X UNFL = DLAMCH( 'Safe minimum' ) X ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' ) X* X* Some Error Checks X* X IF( ITYPE.LT.1 .OR. ITYPE.GT.4 ) THEN X RESULT( 1 ) = TEN / ULP X RETURN X END IF X* X*----------------------------------------------------------------------- X* X* X* Do Test 1 X* X* Norm of A: X* X IF( ITYPE.EQ.4 ) THEN X ANORM = ONE X ELSE X ANORM = MAX( DLANGE( '1', N, N, A, LDA, WORK ), UNFL ) X END IF X* X* Norm of A - UBU' X* X IF( ITYPE.EQ.1 ) THEN X CALL DLACPY( ' ', N, N, A, LDA, WORK, N ) X CALL DGEMM( 'N', 'N', N, N, N, ONE, U, LDU, B, LDB, ZERO, X $ WORK( N**2+1 ), N ) X* X CALL DGEMM( 'N', 'C', N, N, N, -ONE, WORK( N**2+1 ), N, U, LDU, X $ ONE, WORK, N ) X* X ELSE X CALL DLACPY( ' ', N, N, B, LDB, WORK, N ) X* X IF( ITYPE.GE.3 .AND. N.GE.2 ) THEN X CALL DORMC2( 'R', 'T', N, N-1, N-1, U( 2, 1 ), LDU, TAU, X $ WORK( N+1 ), N, WORK( N**2+1 ), IINFO ) X IF( IINFO.NE.0 ) THEN X RESULT( 1 ) = TEN / ULP X RETURN X END IF X* X IF( ITYPE.EQ.3 ) THEN X CALL DORMC2( 'L', 'N', N-1, N, N-1, U( 2, 1 ), LDU, TAU, X $ WORK( 2 ), N, WORK( N**2+1 ), IINFO ) X IF( IINFO.NE.0 ) THEN X RESULT( 1 ) = TEN / ULP X RETURN X END IF X END IF X END IF X* X IF( ITYPE.LT.4 ) THEN X DO 20 JCOL = 1, N X DO 10 JROW = 1, N X WORK( JROW+N*( JCOL-1 ) ) = WORK( JROW+N*( JCOL-1 ) ) X $ - A( JROW, JCOL ) X 10 CONTINUE X 20 CONTINUE X ELSE X DO 30 JDIAG = 1, N X WORK( ( N+1 )*( JDIAG-1 )+1 ) = WORK( ( N+1 )* X $ ( JDIAG-1 )+1 ) - ONE X 30 CONTINUE X END IF X END IF X* X WNORM = DLANGE( '1', N, N, WORK, N, WORK( N**2+1 ) ) X* X IF( ANORM.GT.WNORM ) THEN X RESULT( 1 ) = ( WNORM / ANORM ) / ( N*ULP ) X ELSE X IF( ANORM.LT.ONE ) THEN X RESULT( 1 ) = ( MIN( WNORM, N*ANORM ) / ANORM ) / ( N*ULP ) X ELSE X RESULT( 1 ) = MIN( WNORM / ANORM, DBLE( N ) ) / ( N*ULP ) X END IF X END IF X* X* . . . . . . . . . . . . . . X* X* Do Test 2 X* X* Compute UU' - I X* X IF( ITYPE.EQ.1 ) THEN X CALL DGEMM( 'N', 'C', N, N, N, ONE, U, LDU, U, LDU, ZERO, WORK, X $ N ) X* X DO 40 JDIAG = 1, N X WORK( ( N+1 )*( JDIAG-1 )+1 ) = WORK( ( N+1 )*( JDIAG-1 )+ X $ 1 ) - ONE X 40 CONTINUE X* X RESULT( 2 ) = MIN( DLANGE( '1', N, N, WORK, N, X $ WORK( N**2+1 ) ), DBLE( N ) ) / ( N*ULP ) X END IF X* X*----------------------------------------------------------------------- X* X* X RETURN X* X* End of DGET21 X* X END END_OF_FILE if test 8295 -ne `wc -c <'dget21.f'`; then echo shar: \"'dget21.f'\" unpacked with wrong size! fi # end of 'dget21.f' fi if test -f 'dget22.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'dget22.f'\" else echo shar: Extracting \"'dget22.f'\" \(10597 characters\) sed "s/^X//" >'dget22.f' <<'END_OF_FILE' X SUBROUTINE DGET22( TRANSA, TRANSE, TRANSW, N, A, LDA, E, LDE, WR, X $ WI, WORK, RESULT ) X* X* -- LAPACK test routine (preliminary version) -- X* Univ. of Tennessee, Oak Ridge National Lab, Argonne National Lab, X* Courant Institute, NAG Ltd., and Rice University X* March 26, 1990 X* X* .. Scalar Arguments .. X* X CHARACTER TRANSA, TRANSE, TRANSW X INTEGER LDA, LDE, N X* .. X* X* .. Array Arguments .. X* X DOUBLE PRECISION A( LDA, * ), E( LDE, * ), RESULT( 2 ), WI( * ), X $ WORK( N, * ), WR( * ) X* .. X* X*----------------------------------------------------------------------- X* X* Purpose X* ======= X* X* DGET22 does an eigenvector check. X* X* The basic test is: X* X* RESULT(1) = | A E - E W | / ( |A| |E| ulp ) X* X* using the 1-norm. It also checks the normalization: X* X* RESULT(2) = max | m-norm(E(j)) - 1 | / ( n ulp ) X* j X* X* where E(j) is the j-th eigenvector, and m-norm is the X* max-norm of a vector. If an eigenvector is complex, as X* determined from WI(j) nonzero, then the max-norm of the X* vector ( er + i*ei ) is the maximum of X* |er(1)| + |ei(1)|, ... , |er(n)| + |ei(n)| X* X* W is a block diagonal matrix, with a 1x1 block for each X* real eigenvalue and a 2x2 block for each complex conjugate X* pair. If eigenvalues j and j+1 are a complex conjugate pair, X* so WR(j) = WR(j+1) = wr and WI(j) = - WI(j+1) = wi, then the X* 2 x 2 block corresponding to the pair will be: X* X* ( wr wi ) X* ( -wi wr ) X* X* Such a block multiplying an n x 2 matrix ( ur ui ) on the X* right will be the same as multiplying ur + i*ui by wr + i*wi. X* X* To handle various schemes for storage of left eigenvectors, X* there are options to use (a) A-transpose instead of A, X* (b) E-transpose instead of E, and/or (c) W-transpose instead X* of W. X* X* X* Arguments X* ========== X* X* X* TRANSA - CHARACTER*1 X* If TRANSA='T' or 'C', A-transpose will be used everywhere X* instead of A. If TRANSA='N', A (not transposed) will be X* used. X* Not modified. X* X* TRANSE - CHARACTER*1 X* If TRANSE='T' or 'C', E-transpose will be used everywhere X* instead of E, and the eigenvectors will be in rows of E. X* If TRANSE='N', E (not transposed) will be used, and the X* eigenvectors will be in columns of E. X* Not modified. X* X* TRANSW - CHARACTER*1 X* If TRANSW='T' or 'C', W-transpose will be used everywhere X* instead of W; this corresponds to using -WI(j) instead of X* WI(j) everywhere. If TRANSW='N', W (not transposed) will X* be used. X* Not modified. X* X* N - INTEGER X* The size of the matrix. If it is zero, DGET22 does nothing. X* It must be at least zero. X* Not modified. X* X* A - DOUBLE PRECISION array of dimension ( LDA , N ) X* The matrix whose eigenvectors are in E. X* Not modified. X* X* LDA - INTEGER X* The leading dimension of A. It must be at least 1 X* and at least N. X* Not modified. X* X* E - DOUBLE PRECISION array of dimension ( LDE , N ) X* The matrix of eigenvectors. X* Not modified. X* X* LDE - INTEGER X* The leading dimension of E. It must be at least 1 X* and at least N. X* Not modified. X* X* WR, WI - DOUBLE PRECISION arrays of dimension ( N ). X* The real and imaginary parts of the eigenvalues of A. X* Purely real eigenvalues are indicated by WI(j) = exactly 0. X* Complex conjugate pairs are indicated by WR(j)=WR(j+1) and X* WI(j) = - WI(j+1) non-zero; the real part is assumed to be X* stored in the j-th row/column and the imaginary part in X* the (j+1)-th row/column. These are the only possibilities X* forseen, and strange results may occur if something else X* is supplied. X* Not modified. X* X* WORK - DOUBLE PRECISION array of dimension ( N, N+1 ) X* Workspace. X* Modified. X* X* RESULT - DOUBLE PRECISION array of dimension ( 2 ) X* The value computed by the test described above. X* Modified. X* X*----------------------------------------------------------------------- X* X* X* .. Parameters .. X* X DOUBLE PRECISION ZERO, ONE X PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) X* .. X* X* .. Local Scalars .. X* X CHARACTER NORMA, NORME X INTEGER IECOL, IEROW, INCE, IPAIR, ITRNSE, J, JCOL, X $ JVEC X DOUBLE PRECISION ANORM, ENORM, ENRMAX, ENRMIN, ERRNRM, TEMP1, X $ ULP, UNFL X* .. X* X* .. Local Arrays .. X* X DOUBLE PRECISION WMAT( 2, 2 ) X* .. X* X* .. External Functions .. X* X LOGICAL LSAME X DOUBLE PRECISION DLAMCH, DLANGE X EXTERNAL LSAME, DLAMCH, DLANGE X* .. X* X* .. External Subroutines .. X* X EXTERNAL DAXPY, DGEMM, DLAZRO X* .. X* X* .. Intrinsic Functions .. X* X INTRINSIC ABS, DBLE, MAX, MIN X* .. X* X* X*----------------------------------------------------------------------- X* .. Executable Statements .. X* X* Initialize RESULT (in case N=0) X* X RESULT( 1 ) = ZERO X RESULT( 2 ) = ZERO X IF( N.LE.0 ) X $ RETURN X* X* 1) Constants X* X UNFL = DLAMCH( 'Safe minimum' ) X ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' ) X* X ITRNSE = 0 X INCE = 1 X NORMA = 'O' X NORME = 'O' X* X IF( LSAME( TRANSA, 'T' ) .OR. LSAME( TRANSA, 'C' ) ) THEN X NORMA = 'I' X END IF X IF( LSAME( TRANSE, 'T' ) .OR. LSAME( TRANSE, 'C' ) ) THEN X NORME = 'I' X ITRNSE = 1 X INCE = LDE X END IF X* X*----------------------------------------------------------------------- X* X* Check Normalization of E X* X ENRMIN = ONE / ULP X ENRMAX = ZERO X IF( ITRNSE.EQ.0 ) THEN X* X* . . . . . . . . . . . . . . X* X* Eigenvectors are column vectors. X* X IPAIR = 0 X DO 30 JVEC = 1, N X TEMP1 = ZERO X IF( IPAIR.EQ.0 .AND. JVEC.LT.N .AND. WI( JVEC ).NE.ZERO ) X $ IPAIR = 1 X IF( IPAIR.EQ.1 ) THEN X* X* Complex Eigenvector X* X DO 10 J = 1, N X TEMP1 = MAX( TEMP1, ABS( E( J, JVEC ) )+ X $ ABS( E( J, JVEC+1 ) ) ) X 10 CONTINUE X ENRMIN = MIN( ENRMIN, TEMP1 ) X ENRMAX = MAX( ENRMAX, TEMP1 ) X IPAIR = 2 X ELSE IF( IPAIR.EQ.2 ) THEN X IPAIR = 0 X ELSE X* X* Real Eigenvector X* X DO 20 J = 1, N X TEMP1 = MAX( TEMP1, ABS( E( J, JVEC ) ) ) X 20 CONTINUE X ENRMIN = MIN( ENRMIN, TEMP1 ) X ENRMAX = MAX( ENRMAX, TEMP1 ) X IPAIR = 0 X END IF X 30 CONTINUE X* X* . . . . . . . . . . . . . . X* X* Eigenvectors are row vectors. X* X ELSE X DO 40 JVEC = 1, N X WORK( JVEC, 1 ) = ZERO X 40 CONTINUE X* X DO 60 J = 1, N X IPAIR = 0 X DO 50 JVEC = 1, N X IF( IPAIR.EQ.0 .AND. JVEC.LT.N .AND. WI( JVEC ).NE.ZERO ) X $ IPAIR = 1 X IF( IPAIR.EQ.1 ) THEN X WORK( JVEC, 1 ) = MAX( WORK( JVEC, 1 ), X $ ABS( E( J, JVEC ) )+ X $ ABS( E( J, JVEC+1 ) ) ) X WORK( JVEC+1, 1 ) = WORK( JVEC, 1 ) X ELSE IF( IPAIR.EQ.2 ) THEN X IPAIR = 0 X ELSE X WORK( JVEC, 1 ) = MAX( WORK( JVEC, 1 ), X $ ABS( E( J, JVEC ) ) ) X IPAIR = 0 X END IF X 50 CONTINUE X 60 CONTINUE X* X DO 70 JVEC = 1, N X ENRMIN = MIN( ENRMIN, WORK( JVEC, 1 ) ) X ENRMAX = MAX( ENRMAX, WORK( JVEC, 1 ) ) X 70 CONTINUE X END IF X* X* X* X*----------------------------------------------------------------------- X* X* X* X* Norm of A: X* X ANORM = MAX( DLANGE( NORMA, N, N, A, LDA, WORK ), UNFL ) X* X* Norm of E: X* X ENORM = MAX( DLANGE( NORME, N, N, E, LDE, WORK ), ULP ) X* X* . . . . . . . . . . . . . .. X* X* Norm of Error: X* X* X* Error = AE - EW X* X CALL DLAZRO( N, N, ZERO, ZERO, WORK, N ) X* X IPAIR = 0 X IEROW = 1 X IECOL = 1 X* X DO 80 JCOL = 1, N X IF( ITRNSE.EQ.1 ) THEN X IEROW = JCOL X ELSE X IECOL = JCOL X END IF X* X IF( IPAIR.EQ.0 .AND. WI( JCOL ).NE.ZERO ) X $ IPAIR = 1 X* X IF( IPAIR.EQ.1 ) THEN X WMAT( 1, 1 ) = WR( JCOL ) X WMAT( 2, 1 ) = -WI( JCOL ) X WMAT( 1, 2 ) = WI( JCOL ) X WMAT( 2, 2 ) = WR( JCOL ) X CALL DGEMM( TRANSE, TRANSW, N, 2, 2, ONE, E( IEROW, IECOL ), X $ LDE, WMAT, 2, ZERO, WORK( 1, JCOL ), N ) X IPAIR = 2 X ELSE IF( IPAIR.EQ.2 ) THEN X IPAIR = 0 X* X ELSE X CALL DAXPY( N, WR( JCOL ), E( IEROW, IECOL ), INCE, X $ WORK( 1, JCOL ), 1 ) X IPAIR = 0 X END IF X* X 80 CONTINUE X* X CALL DGEMM( TRANSA, TRANSE, N, N, N, ONE, A, LDA, E, LDE, -ONE, X $ WORK, N ) X* X* X* X ERRNRM = DLANGE( 'One', N, N, WORK, N, WORK( 1, N+1 ) ) / ENORM X* X* . . . . . . . . . . . . . .. X* X* X* Compute RESULT(1) (avoiding under/overflow) X* X* X IF( ANORM.GT.ERRNRM ) THEN X RESULT( 1 ) = ( ERRNRM / ANORM ) / ULP X ELSE X IF( ANORM.LT.ONE ) THEN X RESULT( 1 ) = ( MIN( ERRNRM, ANORM ) / ANORM ) / ULP X ELSE X RESULT( 1 ) = MIN( ERRNRM / ANORM, ONE ) / ULP X END IF X END IF X* X* Compute RESULT(2) : the normalization error in E. X* X RESULT( 2 ) = MAX( ABS( ENRMAX-ONE ), ABS( ENRMIN-ONE ) ) / X $ ( DBLE( N )*ULP ) X* X*----------------------------------------------------------------------- X* X* X RETURN X* X* End of DGET22 X* X END END_OF_FILE if test 10597 -ne `wc -c <'dget22.f'`; then echo shar: \"'dget22.f'\" unpacked with wrong size! fi # end of 'dget22.f' fi if test -f 'dhqr2.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'dhqr2.f'\" else echo shar: Extracting \"'dhqr2.f'\" \(14620 characters\) sed "s/^X//" >'dhqr2.f' <<'END_OF_FILE' X subroutine hqr2(nm,n,low,igh,h,wr,wi,z,ierr) CVD$G noconcur c X integer i,j,k,l,m,n,en,ii,jj,ll,mm,na,nm,nn, X x igh,itn,its,low,mp2,enm2,ierr X double precision h(nm,n),wr(n),wi(n),z(nm,n) X double precision p,q,r,s,t,w,x,y,ra,sa,vi,vr,zz,norm,tst1,tst2 X logical notlas c c this subroutine is a translation of the algol procedure hqr2, c num. math. 16, 181-204(1970) by peters and wilkinson. c handbook for auto. comp., vol.ii-linear algebra, 372-395(1971). c c this subroutine finds the eigenvalues and eigenvectors c of a real upper hessenberg matrix by the qr method. the c eigenvectors of a real general matrix can also be found c if elmhes and eltran or orthes and ortran have c been used to reduce this general matrix to hessenberg form c and to accumulate the similarity transformations. c c on input c c nm must be set to the row dimension of two-dimensional c array parameters as declared in the calling program c dimension statement. c c n is the order of the matrix. c c low and igh are integers determined by the balancing c subroutine balanc. if balanc has not been used, c set low=1, igh=n. c c h contains the upper hessenberg matrix. c c z contains the transformation matrix produced by eltran c after the reduction by elmhes, or by ortran after the c reduction by orthes, if performed. if the eigenvectors c of the hessenberg matrix are desired, z must contain the c identity matrix. c c on output c c h has been destroyed. c c wr and wi contain the real and imaginary parts, c respectively, of the eigenvalues. the eigenvalues c are unordered except that complex conjugate pairs c of values appear consecutively with the eigenvalue c having the positive imaginary part first. if an c error exit is made, the eigenvalues should be correct c for indices ierr+1,...,n. c c z contains the real and imaginary parts of the eigenvectors. c if the i-th eigenvalue is real, the i-th column of z c contains its eigenvector. if the i-th eigenvalue is complex c with positive imaginary part, the i-th and (i+1)-th c columns of z contain the real and imaginary parts of its c eigenvector. the eigenvectors are unnormalized. if an c error exit is made, none of the eigenvectors has been found. c c ierr is set to c zero for normal return, c j if the limit of 30*n iterations is exhausted c while the j-th eigenvalue is being sought. c c calls cdiv for complex division. c c questions and comments should be directed to burton s. garbow, c mathematics and computer science div, argonne national laboratory c c this version dated august 1983. c c ------------------------------------------------------------------ c X ierr = 0 X norm = 0.0d0 X k = 1 c .......... store roots isolated by balanc c and compute matrix norm .......... X do 50 i = 1, n c X do 40 j = k, n X 40 norm = norm + dabs(h(i,j)) c X k = i X if (i .ge. low .and. i .le. igh) go to 50 X wr(i) = h(i,i) X wi(i) = 0.0d0 X 50 continue c X en = igh X t = 0.0d0 X itn = 30*n c .......... search for next eigenvalues .......... X 60 if (en .lt. low) go to 340 X its = 0 X na = en - 1 X enm2 = na - 1 c .......... look for single small sub-diagonal element c for l=en step -1 until low do -- .......... X 70 do 80 ll = low, en X l = en + low - ll X if (l .eq. low) go to 100 X s = dabs(h(l-1,l-1)) + dabs(h(l,l)) X if (s .eq. 0.0d0) s = norm X tst1 = s X tst2 = tst1 + dabs(h(l,l-1)) X if (tst2 .eq. tst1) go to 100 X 80 continue c .......... form shift .......... X 100 x = h(en,en) X if (l .eq. en) go to 270 X y = h(na,na) X w = h(en,na) * h(na,en) X if (l .eq. na) go to 280 X if (itn .eq. 0) go to 1000 X if (its .ne. 10 .and. its .ne. 20) go to 130 c .......... form exceptional shift .......... X t = t + x c X do 120 i = low, en X 120 h(i,i) = h(i,i) - x c X s = dabs(h(en,na)) + dabs(h(na,enm2)) X x = 0.75d0 * s X y = x X w = -0.4375d0 * s * s X 130 its = its + 1 X itn = itn - 1 c .......... look for two consecutive small c sub-diagonal elements. c for m=en-2 step -1 until l do -- .......... X do 140 mm = l, enm2 X m = enm2 + l - mm X zz = h(m,m) X r = x - zz X s = y - zz X p = (r * s - w) / h(m+1,m) + h(m,m+1) X q = h(m+1,m+1) - zz - r - s X r = h(m+2,m+1) X s = dabs(p) + dabs(q) + dabs(r) X p = p / s X q = q / s X r = r / s X if (m .eq. l) go to 150 X tst1 = dabs(p)*(dabs(h(m-1,m-1)) + dabs(zz) + dabs(h(m+1,m+1))) X tst2 = tst1 + dabs(h(m,m-1))*(dabs(q) + dabs(r)) X if (tst2 .eq. tst1) go to 150 X 140 continue c X 150 mp2 = m + 2 c X do 160 i = mp2, en X h(i,i-2) = 0.0d0 X if (i .eq. mp2) go to 160 X h(i,i-3) = 0.0d0 X 160 continue c .......... double qr step involving rows l to en and c columns m to en .......... X do 260 k = m, na X notlas = k .ne. na X if (k .eq. m) go to 170 X p = h(k,k-1) X q = h(k+1,k-1) X r = 0.0d0 X if (notlas) r = h(k+2,k-1) X x = dabs(p) + dabs(q) + dabs(r) X if (x .eq. 0.0d0) go to 260 X p = p / x X q = q / x X r = r / x X 170 s = dsign(dsqrt(p*p+q*q+r*r),p) X if (k .eq. m) go to 180 X h(k,k-1) = -s * x X go to 190 X 180 if (l .ne. m) h(k,k-1) = -h(k,k-1) X 190 p = p + s X x = p / s X y = q / s X zz = r / s X q = q / p X r = r / p X if (notlas) go to 225 c .......... row modification .......... X do 200 j = k, n X p = h(k,j) + q * h(k+1,j) X h(k,j) = h(k,j) - p * x X h(k+1,j) = h(k+1,j) - p * y X 200 continue c X j = min0(en,k+3) c .......... column modification .......... X do 210 i = 1, j X p = x * h(i,k) + y * h(i,k+1) X h(i,k) = h(i,k) - p X h(i,k+1) = h(i,k+1) - p * q X 210 continue c .......... accumulate transformations .......... X do 220 i = low, igh X p = x * z(i,k) + y * z(i,k+1) X z(i,k) = z(i,k) - p X z(i,k+1) = z(i,k+1) - p * q X 220 continue X go to 255 X 225 continue c .......... row modification .......... X do 230 j = k, n X p = h(k,j) + q * h(k+1,j) + r * h(k+2,j) X h(k,j) = h(k,j) - p * x X h(k+1,j) = h(k+1,j) - p * y X h(k+2,j) = h(k+2,j) - p * zz X 230 continue c X j = min0(en,k+3) c .......... column modification .......... X do 240 i = 1, j X p = x * h(i,k) + y * h(i,k+1) + zz * h(i,k+2) X h(i,k) = h(i,k) - p X h(i,k+1) = h(i,k+1) - p * q X h(i,k+2) = h(i,k+2) - p * r X 240 continue c .......... accumulate transformations .......... X do 250 i = low, igh X p = x * z(i,k) + y * z(i,k+1) + zz * z(i,k+2) X z(i,k) = z(i,k) - p X z(i,k+1) = z(i,k+1) - p * q X z(i,k+2) = z(i,k+2) - p * r X 250 continue X 255 continue c X 260 continue c X go to 70 c .......... one root found .......... X 270 h(en,en) = x + t X wr(en) = h(en,en) X wi(en) = 0.0d0 X en = na X go to 60 c .......... two roots found .......... X 280 p = (y - x) / 2.0d0 X q = p * p + w X zz = dsqrt(dabs(q)) X h(en,en) = x + t X x = h(en,en) X h(na,na) = y + t X if (q .lt. 0.0d0) go to 320 c .......... real pair .......... X zz = p + dsign(zz,p) X wr(na) = x + zz X wr(en) = wr(na) X if (zz .ne. 0.0d0) wr(en) = x - w / zz X wi(na) = 0.0d0 X wi(en) = 0.0d0 X x = h(en,na) X s = dabs(x) + dabs(zz) X p = x / s X q = zz / s X r = dsqrt(p*p+q*q) X p = p / r X q = q / r c .......... row modification .......... X do 290 j = na, n X zz = h(na,j) X h(na,j) = q * zz + p * h(en,j) X h(en,j) = q * h(en,j) - p * zz X 290 continue c .......... column modification .......... X do 300 i = 1, en X zz = h(i,na) X h(i,na) = q * zz + p * h(i,en) X h(i,en) = q * h(i,en) - p * zz X 300 continue c .......... accumulate transformations .......... X do 310 i = low, igh X zz = z(i,na) X z(i,na) = q * zz + p * z(i,en) X z(i,en) = q * z(i,en) - p * zz X 310 continue c X go to 330 c .......... complex pair .......... X 320 wr(na) = x + p X wr(en) = x + p X wi(na) = zz X wi(en) = -zz X 330 en = enm2 X go to 60 c .......... all roots found. backsubstitute to find c vectors of upper triangular form .......... X 340 if (norm .eq. 0.0d0) go to 1001 c .......... for en=n step -1 until 1 do -- .......... X do 800 nn = 1, n X en = n + 1 - nn X p = wr(en) X q = wi(en) X na = en - 1 X if (q) 710, 600, 800 c .......... real vector .......... X 600 m = en X h(en,en) = 1.0d0 X if (na .eq. 0) go to 800 c .......... for i=en-1 step -1 until 1 do -- .......... X do 700 ii = 1, na X i = en - ii X w = h(i,i) - p X r = 0.0d0 c X do 610 j = m, en X 610 r = r + h(i,j) * h(j,en) c X if (wi(i) .ge. 0.0d0) go to 630 X zz = w X s = r X go to 700 X 630 m = i X if (wi(i) .ne. 0.0d0) go to 640 X t = w X if (t .ne. 0.0d0) go to 635 X tst1 = norm X t = tst1 X 632 t = 0.01d0 * t X tst2 = norm + t X if (tst2 .gt. tst1) go to 632 X 635 h(i,en) = -r / t X go to 680 c .......... solve real equations .......... X 640 x = h(i,i+1) X y = h(i+1,i) X q = (wr(i) - p) * (wr(i) - p) + wi(i) * wi(i) X t = (x * s - zz * r) / q X h(i,en) = t X if (dabs(x) .le. dabs(zz)) go to 650 X h(i+1,en) = (-r - w * t) / x X go to 680 X 650 h(i+1,en) = (-s - y * t) / zz c c .......... overflow control .......... X 680 t = dabs(h(i,en)) X if (t .eq. 0.0d0) go to 700 X tst1 = t X tst2 = tst1 + 1.0d0/tst1 X if (tst2 .gt. tst1) go to 700 X do 690 j = i, en X h(j,en) = h(j,en)/t X 690 continue c X 700 continue c .......... end real vector .......... X go to 800 c .......... complex vector .......... X 710 m = na c .......... last vector component chosen imaginary so that c eigenvector matrix is triangular .......... X if (dabs(h(en,na)) .le. dabs(h(na,en))) go to 720 X h(na,na) = q / h(en,na) X h(na,en) = -(h(en,en) - p) / h(en,na) X go to 730 X 720 call cdiv(0.0d0,-h(na,en),h(na,na)-p,q,h(na,na),h(na,en)) X 730 h(en,na) = 0.0d0 X h(en,en) = 1.0d0 X enm2 = na - 1 X if (enm2 .eq. 0) go to 800 c .......... for i=en-2 step -1 until 1 do -- .......... X do 795 ii = 1, enm2 X i = na - ii X w = h(i,i) - p X ra = 0.0d0 X sa = 0.0d0 c X do 760 j = m, en X ra = ra + h(i,j) * h(j,na) X sa = sa + h(i,j) * h(j,en) X 760 continue c X if (wi(i) .ge. 0.0d0) go to 770 X zz = w X r = ra X s = sa X go to 795 X 770 m = i X if (wi(i) .ne. 0.0d0) go to 780 X call cdiv(-ra,-sa,w,q,h(i,na),h(i,en)) X go to 790 c .......... solve complex equations .......... X 780 x = h(i,i+1) X y = h(i+1,i) X vr = (wr(i) - p) * (wr(i) - p) + wi(i) * wi(i) - q * q X vi = (wr(i) - p) * 2.0d0 * q X if (vr .ne. 0.0d0 .or. vi .ne. 0.0d0) go to 784 X tst1 = norm * (dabs(w) + dabs(q) + dabs(x) X x + dabs(y) + dabs(zz)) X vr = tst1 X 783 vr = 0.01d0 * vr X tst2 = tst1 + vr X if (tst2 .gt. tst1) go to 783 X 784 call cdiv(x*r-zz*ra+q*sa,x*s-zz*sa-q*ra,vr,vi, X x h(i,na),h(i,en)) X if (dabs(x) .le. dabs(zz) + dabs(q)) go to 785 X h(i+1,na) = (-ra - w * h(i,na) + q * h(i,en)) / x X h(i+1,en) = (-sa - w * h(i,en) - q * h(i,na)) / x X go to 790 X 785 call cdiv(-r-y*h(i,na),-s-y*h(i,en),zz,q, X x h(i+1,na),h(i+1,en)) c c .......... overflow control .......... X 790 t = dmax1(dabs(h(i,na)), dabs(h(i,en))) X if (t .eq. 0.0d0) go to 795 X tst1 = t X tst2 = tst1 + 1.0d0/tst1 X if (tst2 .gt. tst1) go to 795 X do 792 j = i, en X h(j,na) = h(j,na)/t X h(j,en) = h(j,en)/t X 792 continue c X 795 continue c .......... end complex vector .......... X 800 continue c .......... end back substitution. c vectors of isolated roots .......... X do 840 i = 1, n X if (i .ge. low .and. i .le. igh) go to 840 c X do 820 j = i, n X 820 z(i,j) = h(i,j) c X 840 continue c .......... multiply by transformation matrix to give c vectors of original full matrix. c for j=n step -1 until low do -- .......... X do 880 jj = low, n X j = n + low - jj X m = min0(j,igh) c X do 880 i = low, igh X zz = 0.0d0 c X do 860 k = low, m X 860 zz = zz + z(i,k) * h(k,j) c X z(i,j) = zz X 880 continue c X go to 1001 c .......... set error -- all eigenvalues have not c converged after 30*n iterations .......... X 1000 ierr = en X 1001 return X end c X subroutine cdiv(ar,ai,br,bi,cr,ci) X double precision ar,ai,br,bi,cr,ci c c complex division, (cr,ci) = (ar,ai)/(br,bi) c X double precision s,ars,ais,brs,bis X s = dabs(br) + dabs(bi) X ars = ar/s X ais = ai/s X brs = br/s X bis = bi/s X s = brs**2 + bis**2 X cr = (ars*brs + ais*bis)/s X ci = (ais*brs - ars*bis)/s X return X end X END_OF_FILE if test 14620 -ne `wc -c <'dhqr2.f'`; then echo shar: \"'dhqr2.f'\" unpacked with wrong size! fi # end of 'dhqr2.f' fi if test -f 'dhsein.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'dhsein.f'\" else echo shar: Extracting \"'dhsein.f'\" \(16332 characters\) sed "s/^X//" >'dhsein.f' <<'END_OF_FILE' X SUBROUTINE DHSEIN( JOB, SELECT, SOURCE, VECTOR, N, H, LDH, WR, WI, X $ RE, LDRE, LE, LDLE, MM, M, WORK, INFO ) X* X* -- LAPACK routine (preliminary version) -- X* Univ. of Tennessee, Oak Ridge National Lab, Argonne National Lab, X* Courant Institute, NAG Ltd., and Rice University X* March 26, 1990 X* X* .. Scalar Arguments .. X CHARACTER JOB, SOURCE, VECTOR X INTEGER INFO, LDH, LDLE, LDRE, M, MM, N X* .. X* X* .. Array Arguments .. X LOGICAL SELECT( * ) X DOUBLE PRECISION H( LDH, * ), LE( LDLE, * ), RE( LDLE, * ), X $ WI( * ), WORK( * ), WR( * ) X* .. X* X* Purpose X* ======= X* X* This subroutine uses inverse iteration to find specified right X* and/or left eigenvectors of a real upper Hessenberg matrix. X* X* Arguments X* ========= X* X* JOB - CHARACTER*1 X* JOB specifies the computation to be performed by DHSEIN: X* as follows: X* If JOB = 'R', compute right eigenvectors only. X* If JOB = 'L', compute left eigenvectors only. X* If JOB = 'B', compute both right and left eigenvectors. X* Not modified. X* X* SELECT - LOGICAL array, dimension (N). X* SELECT specifies the eigenvectors to be computed. To get X* the eigenvector corresponding to the j-th eigenvalue, set X* SELECT(j) to .TRUE. To get the eigenvectors corresponding X* to a complex conjugate pair of eigenvalues, set the element X* of SELECT corresponding to the first eigenvalue of the pair X* to .TRUE. and the second to .FALSE. (Currently, the value X* of the first element of the pair determines whether the X* pair of eigenvectors is computed.) X* X* On exit, SELECT may have been altered. If the elements of X* SELECT corresponding to a complex conjugate pair of X* eigenvalues were both initially set to .TRUE., the program X* resets the second of the two elements to .FALSE. X* X* X* SOURCE - CHARACTER*1 X* SOURCE specifies the source of eigenvalues supplied in X* WR and WI. X* If SOURCE = 'Q', the eigenvalues were found using DHSEQR; X* thus, if H has zero or negligible sub-diagonal X* entries, and so is block-triangular, then X* the j-th eigenvalue can be assumed to be an X* eigenvalue of the block containing the j-th X* row/column. This property allows DHSEIN to X* perform inverse iteration on just one diagonal X* block. X* If SOURCE = 'N', no assumptions are made on the X* correspondence between eigenvalues and diagonal X* blocks. In this case, DHSEIN must always X* perform inverse iteration using the whole X* matrix H. X* Not modified. X* X* VECTOR - CHARACTER*1 X* VECTOR specifies the source of the initial vectors in X* inverse iteration. X* If VECTOR = 'N', the user does not supply any initial X* vector for the inverse iteration. X* If VECTOR = 'U', the user supplies the initial vectors X* in the array RE and/or LE. The starting vector X* for computing a particular eigenvector will be X* taken from the same place that the eigenvector X* will be stored in. X* Not modified. X* X* N - INTEGER X* The order of the matrix H. X* N must be at least zero. X* Not modified. X* X* H - DOUBLE PRECISION array, dimension (LDH,N). X* H contains the matrix whose eigenvectors are to be computed. X* H must be a real upper Hessenberg matrix. X* Not modified. X* X* LDH - INTEGER. X* The first dimension of H as declared in the calling X* (sub)program. LDH must be at least max(1, N). X* Not modified. X* X* WR,WI - DOUBLE PRECISION arrays, dimension (N). X* On entry, WR and WI contain the real and imaginary parts, X* respectively, of the eigenvalues of H. The N eigenvalues X* may appear in any order, except that a complex conjugate X* pair of eigenvalues must appear consecutively with the X* eigenvalue having the positive imaginary part first. X* Inverse iteration will be performed with each real X* WR(j) for which SELECT(j)=.TRUE. and each complex conjugate X* pair WR(j) +/- i*WI(j) = WR(j+1) -/+ i*WI(j+1) for which X* SELECT(j)=.TRUE. X* X* On exit, WR may have been altered since close eigenvalues X* are perturbed slightly in searching for independent eigen- X* vectors. WI will not be altered. X* X* RE - DOUBLE PRECISION array, dimension (LDRE,MM) X* If VECTOR='U', then on entry RE must contain starting X* vectors for the inverse iteration for the right X* eigenvectors. The starting vector for computing a X* particular eigenvector must be in the same column(s) that X* the eigenvector will be stored in. X* X* The *right* eigenvectors specified by SELECT will be stored X* one after another in the columns of RE, in the same *order* X* (but not necessarily the same position) as their X* eigenvalues. An eigenvector corresponding to a SELECTed X* *real* eigenvalue will take up one column. An eigenvector X* pair corresponding to a SELECTed *complex conjugate pair* X* of eigenvalues will take up two columns: the first column X* will hold the real part, the second will hold the imaginary X* part of the eigenvector corresponding to the eigenvalue X* with *positive* imaginary part. X* X* The eigenvectors will be normalized so that the component X* of largest magnitude is 1; here, the magnitude of a complex X* number x + iy is considered to be |x| + |y|. Eigenvectors X* which do not pass an "acceptance test", i.e., for which the X* inverse iteration does not converge, will be set to zero. X* X* If JOB = 'R' or 'B', RE will be modified. X* If JOB = 'L', RE will not be referenced. X* X* LDRE - INTEGER X* LDRE specifies the leading dimension of RE as declared in X* the calling (sub)program. LDRE must be at least max(1, N). X* If JOB = 'L', LDRE is not referenced. X* Not modified. X* X* LE - DOUBLE PRECISION array, dimension (LDLE,MM) X* If VECTOR='U', then on entry LE must contain starting X* vectors for the inverse iteration for the left eigenvectors. X* The starting vector for computing a particular eigenvector X* must be in the same column(s) that the eigenvector will be X* stored in. X* X* The conjugate transposes of the *left* eigenvectors X* specified by SELECT will be stored one after another in the X* columns of LE, in the same *order* (but not necessarily the X* same position) as their eigenvalues. An eigenvector X* corresponding to a SELECTed *real* eigenvalue will take up X* one column. An eigenvector pair corresponding to a X* SELECTed *complex conjugate pair* of eigenvalues will take X* up two columns: the first column will hold the real part, X* the second will hold the imaginary part of the conjugate X* transpose of the left eigenvector corresponding to the X* eigenvalue with *positive* imaginary part. X* X* The eigenvectors will be normalized so that the component X* of largest magnitude is 1; here, the magnitude of a complex X* number x + iy is considered to be |x| + |y|. Eigenvectors X* which do not pass an "acceptance test", i.e., for which the X* inverse iteration does not converge, will be set to zero. X* X* If JOB = 'L' or 'B', LE will be modified. X* If JOB = 'R', LE will not be referenced. X* X* X* LDLE - INTEGER X* LDLE specifies the leading dimension of LE as declared in X* the calling (sub)program. LDLE must be at least max(1, N). X* If JOB = 'R', LDLE is not referenced. X* Not modified. X* X* MM - INTEGER X* The number of columns in LE and/or RE. Note that X* two columns are required to store the eigenvector X* corresponding to a complex eigenvalue. X* Not modified. X* X* M - INTEGER X* On exit, M is the number of columns in LE and/or RE actually X* used to store the eigenvectors. X* X* WORK - DOUBLE PRECISION array, dimension ( (N+2)**2 + N ). X* WORK is a (N+2)**2 + N workarray. X* X* INFO - INTEGER X* On exit, INFO is set to X* 0 for normal return, X* -k if input argument number k is illegal. X* N+1 if more than MM columns of RE and/or LE are X* necessary to store the SELECTed eigenvectors. X* X* .. Parameters .. X DOUBLE PRECISION ZERO, ONE X PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) X* .. X* X* .. Local Scalars .. X INTEGER I, IJOB, IP, IP1, ISOURC, IVECTO, K, LK, LK1, X $ S, UK, UK1 X DOUBLE PRECISION BIGNUM, EPS3, ILAMBD, NORM, NORMF, NORML, X $ NORMR, OVFL, RLAMBD, SMLNUM, ULP, UNFL X* .. X* X* .. External Functions .. X LOGICAL LSAME X DOUBLE PRECISION DLAMCH, DLANHS X EXTERNAL LSAME, DLAMCH, DLANHS X* .. X* X* .. External Subroutines .. X EXTERNAL DLAEIN, XERBLA X* .. X* X* .. Intrinsic Functions .. X INTRINSIC ABS, MAX X* .. X* X* .. Executable Statements .. X* X* Deconde and Test the input parameter X* X IF( LSAME( JOB, 'R' ) ) THEN X IJOB = 1 X ELSE IF( LSAME( JOB, 'L' ) ) THEN X IJOB = 2 X ELSE IF( LSAME( JOB, 'B' ) ) THEN X IJOB = 3 X ELSE X IJOB = -1 X END IF X* X IF( LSAME( SOURCE, 'Q' ) ) THEN X ISOURC = 1 X ELSE IF( LSAME( SOURCE, 'N' ) ) THEN X ISOURC = 2 X ELSE X ISOURC = -1 X END IF X* X IF( LSAME( VECTOR, 'N' ) ) THEN X IVECTO = 1 X ELSE IF( LSAME( VECTOR, 'U' ) ) THEN X IVECTO = 2 X ELSE X IVECTO = -1 X END IF X* X INFO = 0 X IF( IJOB.EQ.-1 ) THEN X INFO = -1 X ELSE IF( ISOURC.EQ.-1 ) THEN X INFO = -3 X ELSE IF( IVECTO.EQ.-1 ) THEN X INFO = -4 X ELSE IF( N.LT.0 ) THEN X INFO = -5 X ELSE IF( LDH.LT.MAX( 1, N ) ) THEN X INFO = -7 X END IF X IF( IJOB.EQ.1 .OR. IJOB.EQ.3 ) THEN X IF( LDRE.LT.MAX( 1, N ) ) X $ INFO = -11 X END IF X IF( IJOB.EQ.2 .OR. IJOB.EQ.3 ) THEN X IF( LDLE.LT.MAX( 1, N ) ) X $ INFO = -13 X END IF X IF( INFO.NE.0 ) THEN X CALL XERBLA( 'DHSEIN', -INFO ) X RETURN X END IF X* X* Quick return if possible X* X IF( N.EQ.0 ) X $ RETURN X* X* Set constants to control overflow. X* X UNFL = DLAMCH( 'Safe minimum' ) X OVFL = DLAMCH( 'Overflow' ) X ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' ) X SMLNUM = MAX( UNFL*( N/ULP ), N/( ULP*OVFL ) ) X BIGNUM = ( ONE-ULP ) / SMLNUM X* X* Compute inf-norm of matrix H. X* X NORMF = DLANHS( 'I', N, H, LDH, WORK ) X* X* ip = 0, real eigenvalue, X* 1, first of conjugate complex pair, X* -1, second of conjugate complex pair. X* X UK = 0 X LK = 1 X IP = 0 X S = 1 X* X DO 170 K = 1, N X IF( IP.EQ.-1 ) X $ GO TO 160 X IF( WI( K ).EQ.ZERO ) X $ GO TO 10 X IP = 1 X IF( SELECT( K ) .AND. SELECT( K+1 ) ) X $ SELECT( K+1 ) = .FALSE. X 10 CONTINUE X IF( .NOT.SELECT( K ) ) X $ GO TO 160 X IF( IP.EQ.0 .AND. S.GT.MM ) X $ GO TO 180 X IF( IP.NE.0 .AND. S+1.GT.MM ) X $ GO TO 180 X* X* If the affiliation of eigenvalue is known, split checking X* X IF( ISOURC.EQ.1 ) THEN X* X IF( IJOB.EQ.1 .OR. IJOB.EQ.3 ) THEN X* X* Split checking for right eigenvector. The inverse X* iteration works on H(1:UK,1:UK) X* X IF( UK.GE.K ) X $ GO TO 50 X DO 20 UK1 = K, N X IF( UK1.EQ.N ) X $ GO TO 40 X IF( H( UK1+1, UK1 ).EQ.ZERO ) X $ GO TO 30 X 20 CONTINUE X* X 30 CONTINUE X UK = UK1 X NORMR = DLANHS( 'I', UK, H, LDH, WORK ) X GO TO 50 X 40 CONTINUE X UK = N X NORMR = NORMF X END IF X* X* Split checking for left eigenvector. The inverse X* iteration works on H(LK:N,LK:N). X* X 50 CONTINUE X IF( IJOB.EQ.2 .OR. IJOB.EQ.3 ) THEN X DO 60 LK1 = K, LK, -1 X IF( LK1.EQ.LK ) X $ GO TO 70 X IF( H( LK1, LK1-1 ).EQ.ZERO ) X $ GO TO 70 X 60 CONTINUE X* X 70 CONTINUE X IF( LK1.EQ.1 ) THEN X LK = LK1 X NORML = NORMF X GO TO 80 X END IF X* X IF( LK1.EQ.LK ) THEN X LK = LK1 X GO TO 80 X ELSE X LK = LK1 X NORML = DLANHS( 'I', N-LK+1, H( LK, LK ), LDH, WORK ) X END IF X END IF X* X 80 CONTINUE X IF( IJOB.EQ.1 ) THEN X NORM = NORMR X ELSE IF( IJOB.EQ.2 ) THEN X NORM = NORML X ELSE IF( IJOB.EQ.3 ) THEN X NORM = MAX( NORMR, NORML ) X END IF X* X ELSE X* X* If the affiliation of eigenvalue is not known, the X*