octave: libcruft/lapack/claqr0.f annotate

author	David Bateman <dbateman@free.fr>
date	Sun, 27 Apr 2008 22:34:17 +0200
parents
children

rev	line source
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	1 SUBROUTINE CLAQR0( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	2 $ IHIZ, Z, LDZ, WORK, LWORK, INFO )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	3 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	4 * -- LAPACK auxiliary routine (version 3.1) --
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	6 * November 2006
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	7 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	8 * .. Scalar Arguments ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	9 INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, LWORK, N
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	10 LOGICAL WANTT, WANTZ
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	11 * ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	12 * .. Array Arguments ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	13 COMPLEX H( LDH, * ), W( * ), WORK( * ), Z( LDZ, * )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	14 * ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	15 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	16 * Purpose
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	17 * =======
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	18 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	19 * CLAQR0 computes the eigenvalues of a Hessenberg matrix H
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	20 * and, optionally, the matrices T and Z from the Schur decomposition
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	21 * H = Z T Z**H, where T is an upper triangular matrix (the
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	22 * Schur form), and Z is the unitary matrix of Schur vectors.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	23 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	24 * Optionally Z may be postmultiplied into an input unitary
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	25 * matrix Q so that this routine can give the Schur factorization
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	26 * of a matrix A which has been reduced to the Hessenberg form H
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	27 * by the unitary matrix Q: A = QHQ*H = (QZ)H(QZ)*H.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	28 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	29 * Arguments
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	30 * =========
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	31 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	32 * WANTT (input) LOGICAL
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	33 * = .TRUE. : the full Schur form T is required;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	34 * = .FALSE.: only eigenvalues are required.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	35 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	36 * WANTZ (input) LOGICAL
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	37 * = .TRUE. : the matrix of Schur vectors Z is required;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	38 * = .FALSE.: Schur vectors are not required.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	39 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	40 * N (input) INTEGER
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	41 * The order of the matrix H. N .GE. 0.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	42 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	43 * ILO (input) INTEGER
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	44 * IHI (input) INTEGER
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	45 * It is assumed that H is already upper triangular in rows
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	46 * and columns 1:ILO-1 and IHI+1:N and, if ILO.GT.1,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	47 * H(ILO,ILO-1) is zero. ILO and IHI are normally set by a
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	48 * previous call to CGEBAL, and then passed to CGEHRD when the
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	49 * matrix output by CGEBAL is reduced to Hessenberg form.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	50 * Otherwise, ILO and IHI should be set to 1 and N,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	51 * respectively. If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	52 * If N = 0, then ILO = 1 and IHI = 0.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	53 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	54 * H (input/output) COMPLEX array, dimension (LDH,N)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	55 * On entry, the upper Hessenberg matrix H.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	56 * On exit, if INFO = 0 and WANTT is .TRUE., then H
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	57 * contains the upper triangular matrix T from the Schur
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	58 * decomposition (the Schur form). If INFO = 0 and WANT is
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	59 * .FALSE., then the contents of H are unspecified on exit.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	60 * (The output value of H when INFO.GT.0 is given under the
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	61 * description of INFO below.)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	62 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	63 * This subroutine may explicitly set H(i,j) = 0 for i.GT.j and
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	64 * j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	65 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	66 * LDH (input) INTEGER
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	67 * The leading dimension of the array H. LDH .GE. max(1,N).
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	68 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	69 * W (output) COMPLEX array, dimension (N)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	70 * The computed eigenvalues of H(ILO:IHI,ILO:IHI) are stored
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	71 * in W(ILO:IHI). If WANTT is .TRUE., then the eigenvalues are
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	72 * stored in the same order as on the diagonal of the Schur
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	73 * form returned in H, with W(i) = H(i,i).
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	74 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	75 * Z (input/output) COMPLEX array, dimension (LDZ,IHI)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	76 * If WANTZ is .FALSE., then Z is not referenced.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	77 * If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	78 * replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	79 * orthogonal Schur factor of H(ILO:IHI,ILO:IHI).
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	80 * (The output value of Z when INFO.GT.0 is given under
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	81 * the description of INFO below.)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	82 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	83 * LDZ (input) INTEGER
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	84 * The leading dimension of the array Z. if WANTZ is .TRUE.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	85 * then LDZ.GE.MAX(1,IHIZ). Otherwize, LDZ.GE.1.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	86 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	87 * WORK (workspace/output) COMPLEX array, dimension LWORK
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	88 * On exit, if LWORK = -1, WORK(1) returns an estimate of
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	89 * the optimal value for LWORK.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	90 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	91 * LWORK (input) INTEGER
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	92 * The dimension of the array WORK. LWORK .GE. max(1,N)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	93 * is sufficient, but LWORK typically as large as 6*N may
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	94 * be required for optimal performance. A workspace query
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	95 * to determine the optimal workspace size is recommended.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	96 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	97 * If LWORK = -1, then CLAQR0 does a workspace query.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	98 * In this case, CLAQR0 checks the input parameters and
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	99 * estimates the optimal workspace size for the given
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	100 * values of N, ILO and IHI. The estimate is returned
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	101 * in WORK(1). No error message related to LWORK is
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	102 * issued by XERBLA. Neither H nor Z are accessed.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	103 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	104 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	105 * INFO (output) INTEGER
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	106 * = 0: successful exit
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	107 * .GT. 0: if INFO = i, CLAQR0 failed to compute all of
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	108 * the eigenvalues. Elements 1:ilo-1 and i+1:n of WR
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	109 * and WI contain those eigenvalues which have been
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	110 * successfully computed. (Failures are rare.)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	111 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	112 * If INFO .GT. 0 and WANT is .FALSE., then on exit,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	113 * the remaining unconverged eigenvalues are the eigen-
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	114 * values of the upper Hessenberg matrix rows and
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	115 * columns ILO through INFO of the final, output
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	116 * value of H.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	117 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	118 * If INFO .GT. 0 and WANTT is .TRUE., then on exit
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	119 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	120 * () (initial value of H)U = U*(final value of H)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	121 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	122 * where U is a unitary matrix. The final
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	123 * value of H is upper Hessenberg and triangular in
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	124 * rows and columns INFO+1 through IHI.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	125 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	126 * If INFO .GT. 0 and WANTZ is .TRUE., then on exit
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	127 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	128 * (final value of Z(ILO:IHI,ILOZ:IHIZ)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	129 * = (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	130 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	131 * where U is the unitary matrix in (*) (regard-
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	132 * less of the value of WANTT.)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	133 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	134 * If INFO .GT. 0 and WANTZ is .FALSE., then Z is not
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	135 * accessed.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	136 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	137 * ================================================================
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	138 * Based on contributions by
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	139 * Karen Braman and Ralph Byers, Department of Mathematics,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	140 * University of Kansas, USA
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	141 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	142 * ================================================================
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	143 * References:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	144 * K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	145 * Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	146 * Performance, SIAM Journal of Matrix Analysis, volume 23, pages
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	147 * 929--947, 2002.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	148 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	149 * K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	150 * Algorithm Part II: Aggressive Early Deflation, SIAM Journal
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	151 * of Matrix Analysis, volume 23, pages 948--973, 2002.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	152 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	153 * ================================================================
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	154 * .. Parameters ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	155 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	156 * ==== Matrices of order NTINY or smaller must be processed by
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	157 * . CLAHQR because of insufficient subdiagonal scratch space.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	158 * . (This is a hard limit.) ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	159 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	160 * ==== Exceptional deflation windows: try to cure rare
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	161 * . slow convergence by increasing the size of the
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	162 * . deflation window after KEXNW iterations. =====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	163 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	164 * ==== Exceptional shifts: try to cure rare slow convergence
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	165 * . with ad-hoc exceptional shifts every KEXSH iterations.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	166 * . The constants WILK1 and WILK2 are used to form the
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	167 * . exceptional shifts. ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	168 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	169 INTEGER NTINY
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	170 PARAMETER ( NTINY = 11 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	171 INTEGER KEXNW, KEXSH
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	172 PARAMETER ( KEXNW = 5, KEXSH = 6 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	173 REAL WILK1
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	174 PARAMETER ( WILK1 = 0.75e0 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	175 COMPLEX ZERO, ONE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	176 PARAMETER ( ZERO = ( 0.0e0, 0.0e0 ),
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	177 $ ONE = ( 1.0e0, 0.0e0 ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	178 REAL TWO
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	179 PARAMETER ( TWO = 2.0e0 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	180 * ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	181 * .. Local Scalars ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	182 COMPLEX AA, BB, CC, CDUM, DD, DET, RTDISC, SWAP, TR2
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	183 REAL S
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	184 INTEGER I, INF, IT, ITMAX, K, KACC22, KBOT, KDU, KS,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	185 $ KT, KTOP, KU, KV, KWH, KWTOP, KWV, LD, LS,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	186 $ LWKOPT, NDFL, NH, NHO, NIBBLE, NMIN, NS, NSMAX,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	187 $ NSR, NVE, NW, NWMAX, NWR
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	188 LOGICAL NWINC, SORTED
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	189 CHARACTER JBCMPZ*2
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	190 * ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	191 * .. External Functions ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	192 INTEGER ILAENV
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	193 EXTERNAL ILAENV
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	194 * ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	195 * .. Local Arrays ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	196 COMPLEX ZDUM( 1, 1 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	197 * ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	198 * .. External Subroutines ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	199 EXTERNAL CLACPY, CLAHQR, CLAQR3, CLAQR4, CLAQR5
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	200 * ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	201 * .. Intrinsic Functions ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	202 INTRINSIC ABS, AIMAG, CMPLX, INT, MAX, MIN, MOD, REAL,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	203 $ SQRT
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	204 * ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	205 * .. Statement Functions ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	206 REAL CABS1
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	207 * ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	208 * .. Statement Function definitions ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	209 CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	210 * ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	211 * .. Executable Statements ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	212 INFO = 0
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	213 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	214 * ==== Quick return for N = 0: nothing to do. ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	215 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	216 IF( N.EQ.0 ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	217 WORK( 1 ) = ONE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	218 RETURN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	219 END IF
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	220 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	221 * ==== Set up job flags for ILAENV. ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	222 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	223 IF( WANTT ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	224 JBCMPZ( 1: 1 ) = 'S'
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	225 ELSE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	226 JBCMPZ( 1: 1 ) = 'E'
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	227 END IF
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	228 IF( WANTZ ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	229 JBCMPZ( 2: 2 ) = 'V'
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	230 ELSE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	231 JBCMPZ( 2: 2 ) = 'N'
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	232 END IF
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	233 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	234 * ==== Tiny matrices must use CLAHQR. ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	235 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	236 IF( N.LE.NTINY ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	237 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	238 * ==== Estimate optimal workspace. ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	239 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	240 LWKOPT = 1
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	241 IF( LWORK.NE.-1 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	242 $ CALL CLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	243 $ IHIZ, Z, LDZ, INFO )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	244 ELSE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	245 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	246 * ==== Use small bulge multi-shift QR with aggressive early
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	247 * . deflation on larger-than-tiny matrices. ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	248 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	249 * ==== Hope for the best. ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	250 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	251 INFO = 0
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	252 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	253 * ==== NWR = recommended deflation window size. At this
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	254 * . point, N .GT. NTINY = 11, so there is enough
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	255 * . subdiagonal workspace for NWR.GE.2 as required.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	256 * . (In fact, there is enough subdiagonal space for
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	257 * . NWR.GE.3.) ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	258 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	259 NWR = ILAENV( 13, 'CLAQR0', JBCMPZ, N, ILO, IHI, LWORK )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	260 NWR = MAX( 2, NWR )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	261 NWR = MIN( IHI-ILO+1, ( N-1 ) / 3, NWR )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	262 NW = NWR
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	263 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	264 * ==== NSR = recommended number of simultaneous shifts.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	265 * . At this point N .GT. NTINY = 11, so there is at
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	266 * . enough subdiagonal workspace for NSR to be even
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	267 * . and greater than or equal to two as required. ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	268 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	269 NSR = ILAENV( 15, 'CLAQR0', JBCMPZ, N, ILO, IHI, LWORK )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	270 NSR = MIN( NSR, ( N+6 ) / 9, IHI-ILO )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	271 NSR = MAX( 2, NSR-MOD( NSR, 2 ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	272 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	273 * ==== Estimate optimal workspace ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	274 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	275 * ==== Workspace query call to CLAQR3 ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	276 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	277 CALL CLAQR3( WANTT, WANTZ, N, ILO, IHI, NWR+1, H, LDH, ILOZ,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	278 $ IHIZ, Z, LDZ, LS, LD, W, H, LDH, N, H, LDH, N, H,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	279 $ LDH, WORK, -1 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	280 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	281 * ==== Optimal workspace = MAX(CLAQR5, CLAQR3) ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	282 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	283 LWKOPT = MAX( 3*NSR / 2, INT( WORK( 1 ) ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	284 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	285 * ==== Quick return in case of workspace query. ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	286 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	287 IF( LWORK.EQ.-1 ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	288 WORK( 1 ) = CMPLX( LWKOPT, 0 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	289 RETURN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	290 END IF
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	291 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	292 * ==== CLAHQR/CLAQR0 crossover point ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	293 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	294 NMIN = ILAENV( 12, 'CLAQR0', JBCMPZ, N, ILO, IHI, LWORK )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	295 NMIN = MAX( NTINY, NMIN )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	296 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	297 * ==== Nibble crossover point ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	298 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	299 NIBBLE = ILAENV( 14, 'CLAQR0', JBCMPZ, N, ILO, IHI, LWORK )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	300 NIBBLE = MAX( 0, NIBBLE )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	301 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	302 * ==== Accumulate reflections during ttswp? Use block
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	303 * . 2-by-2 structure during matrix-matrix multiply? ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	304 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	305 KACC22 = ILAENV( 16, 'CLAQR0', JBCMPZ, N, ILO, IHI, LWORK )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	306 KACC22 = MAX( 0, KACC22 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	307 KACC22 = MIN( 2, KACC22 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	308 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	309 * ==== NWMAX = the largest possible deflation window for
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	310 * . which there is sufficient workspace. ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	311 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	312 NWMAX = MIN( ( N-1 ) / 3, LWORK / 2 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	313 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	314 * ==== NSMAX = the Largest number of simultaneous shifts
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	315 * . for which there is sufficient workspace. ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	316 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	317 NSMAX = MIN( ( N+6 ) / 9, 2*LWORK / 3 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	318 NSMAX = NSMAX - MOD( NSMAX, 2 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	319 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	320 * ==== NDFL: an iteration count restarted at deflation. ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	321 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	322 NDFL = 1
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	323 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	324 * ==== ITMAX = iteration limit ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	325 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	326 ITMAX = MAX( 30, 2KEXSH )MAX( 10, ( IHI-ILO+1 ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	327 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	328 * ==== Last row and column in the active block ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	329 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	330 KBOT = IHI
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	331 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	332 * ==== Main Loop ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	333 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	334 DO 70 IT = 1, ITMAX
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	335 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	336 * ==== Done when KBOT falls below ILO ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	337 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	338 IF( KBOT.LT.ILO )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	339 $ GO TO 80
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	340 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	341 * ==== Locate active block ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	342 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	343 DO 10 K = KBOT, ILO + 1, -1
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	344 IF( H( K, K-1 ).EQ.ZERO )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	345 $ GO TO 20
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	346 10 CONTINUE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	347 K = ILO
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	348 20 CONTINUE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	349 KTOP = K
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	350 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	351 * ==== Select deflation window size ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	352 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	353 NH = KBOT - KTOP + 1
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	354 IF( NDFL.LT.KEXNW .OR. NH.LT.NW ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	355 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	356 * ==== Typical deflation window. If possible and
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	357 * . advisable, nibble the entire active block.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	358 * . If not, use size NWR or NWR+1 depending upon
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	359 * . which has the smaller corresponding subdiagonal
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	360 * . entry (a heuristic). ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	361 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	362 NWINC = .TRUE.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	363 IF( NH.LE.MIN( NMIN, NWMAX ) ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	364 NW = NH
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	365 ELSE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	366 NW = MIN( NWR, NH, NWMAX )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	367 IF( NW.LT.NWMAX ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	368 IF( NW.GE.NH-1 ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	369 NW = NH
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	370 ELSE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	371 KWTOP = KBOT - NW + 1
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	372 IF( CABS1( H( KWTOP, KWTOP-1 ) ).GT.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	373 $ CABS1( H( KWTOP-1, KWTOP-2 ) ) )NW = NW + 1
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	374 END IF
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	375 END IF
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	376 END IF
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	377 ELSE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	378 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	379 * ==== Exceptional deflation window. If there have
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	380 * . been no deflations in KEXNW or more iterations,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	381 * . then vary the deflation window size. At first,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	382 * . because, larger windows are, in general, more
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	383 * . powerful than smaller ones, rapidly increase the
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	384 * . window up to the maximum reasonable and possible.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	385 * . Then maybe try a slightly smaller window. ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	386 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	387 IF( NWINC .AND. NW.LT.MIN( NWMAX, NH ) ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	388 NW = MIN( NWMAX, NH, 2*NW )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	389 ELSE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	390 NWINC = .FALSE.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	391 IF( NW.EQ.NH .AND. NH.GT.2 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	392 $ NW = NH - 1
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	393 END IF
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	394 END IF
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	395 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	396 * ==== Aggressive early deflation:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	397 * . split workspace under the subdiagonal into
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	398 * . - an nw-by-nw work array V in the lower
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	399 * . left-hand-corner,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	400 * . - an NW-by-at-least-NW-but-more-is-better
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	401 * . (NW-by-NHO) horizontal work array along
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	402 * . the bottom edge,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	403 * . - an at-least-NW-but-more-is-better (NHV-by-NW)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	404 * . vertical work array along the left-hand-edge.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	405 * . ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	406 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	407 KV = N - NW + 1
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	408 KT = NW + 1
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	409 NHO = ( N-NW-1 ) - KT + 1
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	410 KWV = NW + 2
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	411 NVE = ( N-NW ) - KWV + 1
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	412 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	413 * ==== Aggressive early deflation ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	414 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	415 CALL CLAQR3( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	416 $ IHIZ, Z, LDZ, LS, LD, W, H( KV, 1 ), LDH, NHO,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	417 $ H( KV, KT ), LDH, NVE, H( KWV, 1 ), LDH, WORK,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	418 $ LWORK )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	419 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	420 * ==== Adjust KBOT accounting for new deflations. ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	421 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	422 KBOT = KBOT - LD
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	423 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	424 * ==== KS points to the shifts. ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	425 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	426 KS = KBOT - LS + 1
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	427 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	428 * ==== Skip an expensive QR sweep if there is a (partly
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	429 * . heuristic) reason to expect that many eigenvalues
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	430 * . will deflate without it. Here, the QR sweep is
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	431 * . skipped if many eigenvalues have just been deflated
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	432 * . or if the remaining active block is small.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	433 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	434 IF( ( LD.EQ.0 ) .OR. ( ( 100LD.LE.NWNIBBLE ) .AND. ( KBOT-
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	435 $ KTOP+1.GT.MIN( NMIN, NWMAX ) ) ) ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	436 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	437 * ==== NS = nominal number of simultaneous shifts.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	438 * . This may be lowered (slightly) if CLAQR3
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	439 * . did not provide that many shifts. ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	440 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	441 NS = MIN( NSMAX, NSR, MAX( 2, KBOT-KTOP ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	442 NS = NS - MOD( NS, 2 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	443 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	444 * ==== If there have been no deflations
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	445 * . in a multiple of KEXSH iterations,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	446 * . then try exceptional shifts.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	447 * . Otherwise use shifts provided by
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	448 * . CLAQR3 above or from the eigenvalues
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	449 * . of a trailing principal submatrix. ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	450 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	451 IF( MOD( NDFL, KEXSH ).EQ.0 ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	452 KS = KBOT - NS + 1
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	453 DO 30 I = KBOT, KS + 1, -2
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	454 W( I ) = H( I, I ) + WILK1*CABS1( H( I, I-1 ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	455 W( I-1 ) = W( I )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	456 30 CONTINUE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	457 ELSE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	458 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	459 * ==== Got NS/2 or fewer shifts? Use CLAQR4 or
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	460 * . CLAHQR on a trailing principal submatrix to
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	461 * . get more. (Since NS.LE.NSMAX.LE.(N+6)/9,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	462 * . there is enough space below the subdiagonal
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	463 * . to fit an NS-by-NS scratch array.) ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	464 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	465 IF( KBOT-KS+1.LE.NS / 2 ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	466 KS = KBOT - NS + 1
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	467 KT = N - NS + 1
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	468 CALL CLACPY( 'A', NS, NS, H( KS, KS ), LDH,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	469 $ H( KT, 1 ), LDH )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	470 IF( NS.GT.NMIN ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	471 CALL CLAQR4( .false., .false., NS, 1, NS,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	472 $ H( KT, 1 ), LDH, W( KS ), 1, 1,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	473 $ ZDUM, 1, WORK, LWORK, INF )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	474 ELSE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	475 CALL CLAHQR( .false., .false., NS, 1, NS,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	476 $ H( KT, 1 ), LDH, W( KS ), 1, 1,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	477 $ ZDUM, 1, INF )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	478 END IF
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	479 KS = KS + INF
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	480 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	481 * ==== In case of a rare QR failure use
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	482 * . eigenvalues of the trailing 2-by-2
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	483 * . principal submatrix. Scale to avoid
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	484 * . overflows, underflows and subnormals.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	485 * . (The scale factor S can not be zero,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	486 * . because H(KBOT,KBOT-1) is nonzero.) ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	487 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	488 IF( KS.GE.KBOT ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	489 S = CABS1( H( KBOT-1, KBOT-1 ) ) +
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	490 $ CABS1( H( KBOT, KBOT-1 ) ) +
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	491 $ CABS1( H( KBOT-1, KBOT ) ) +
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	492 $ CABS1( H( KBOT, KBOT ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	493 AA = H( KBOT-1, KBOT-1 ) / S
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	494 CC = H( KBOT, KBOT-1 ) / S
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	495 BB = H( KBOT-1, KBOT ) / S
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	496 DD = H( KBOT, KBOT ) / S
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	497 TR2 = ( AA+DD ) / TWO
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	498 DET = ( AA-TR2 )( DD-TR2 ) - BBCC
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	499 RTDISC = SQRT( -DET )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	500 W( KBOT-1 ) = ( TR2+RTDISC )*S

7789

82be108cc558 First attempt at single precision tyeps