octave-nkf: libcruft/lapack/claqr5.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 CLAQR5( WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS, S,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	2 $ H, LDH, ILOZ, IHIZ, Z, LDZ, V, LDV, U, LDU, NV,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	3 $ WV, LDWV, NH, WH, LDWH )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	4 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	5 * -- LAPACK auxiliary routine (version 3.1) --
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	6 * 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	7 * November 2006
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	8 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	9 * .. Scalar Arguments ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	10 INTEGER IHIZ, ILOZ, KACC22, KBOT, KTOP, LDH, LDU, LDV,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	11 $ LDWH, LDWV, LDZ, N, NH, NSHFTS, NV
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	12 LOGICAL WANTT, WANTZ
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	13 * ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	14 * .. Array Arguments ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	15 COMPLEX H( LDH, * ), S( * ), U( LDU, * ), V( LDV, * ),
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	16 $ WH( LDWH, * ), WV( LDWV, * ), Z( LDZ, * )
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 * This auxiliary subroutine called by CLAQR0 performs a
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	20 * single small-bulge multi-shift QR sweep.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	21 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	22 * WANTT (input) logical scalar
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	23 * WANTT = .true. if the triangular Schur factor
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	24 * is being computed. WANTT is set to .false. otherwise.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	25 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	26 * WANTZ (input) logical scalar
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	27 * WANTZ = .true. if the unitary Schur factor is being
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	28 * computed. WANTZ is set to .false. otherwise.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	29 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	30 * KACC22 (input) integer with value 0, 1, or 2.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	31 * Specifies the computation mode of far-from-diagonal
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	32 * orthogonal updates.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	33 * = 0: CLAQR5 does not accumulate reflections and does not
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	34 * use matrix-matrix multiply to update far-from-diagonal
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	35 * matrix entries.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	36 * = 1: CLAQR5 accumulates reflections and uses matrix-matrix
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	37 * multiply to update the far-from-diagonal matrix entries.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	38 * = 2: CLAQR5 accumulates reflections, uses matrix-matrix
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	39 * multiply to update the far-from-diagonal matrix entries,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	40 * and takes advantage of 2-by-2 block structure during
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	41 * matrix multiplies.
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 * N (input) integer scalar
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	44 * N is the order of the Hessenberg matrix H upon which this
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	45 * subroutine operates.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	46 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	47 * KTOP (input) integer scalar
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	48 * KBOT (input) integer scalar
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	49 * These are the first and last rows and columns of an
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	50 * isolated diagonal block upon which the QR sweep is to be
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	51 * applied. It is assumed without a check that
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	52 * either KTOP = 1 or H(KTOP,KTOP-1) = 0
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	53 * and
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	54 * either KBOT = N or H(KBOT+1,KBOT) = 0.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	55 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	56 * NSHFTS (input) integer scalar
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	57 * NSHFTS gives the number of simultaneous shifts. NSHFTS
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	58 * must be positive and even.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	59 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	60 * S (input) COMPLEX array of size (NSHFTS)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	61 * S contains the shifts of origin that define the multi-
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	62 * shift QR sweep.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	63 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	64 * H (input/output) COMPLEX array of size (LDH,N)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	65 * On input H contains a Hessenberg matrix. On output a
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	66 * multi-shift QR sweep with shifts SR(J)+i*SI(J) is applied
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	67 * to the isolated diagonal block in rows and columns KTOP
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	68 * through KBOT.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	69 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	70 * LDH (input) integer scalar
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	71 * LDH is the leading dimension of H just as declared in the
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	72 * calling procedure. LDH.GE.MAX(1,N).
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	73 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	74 * ILOZ (input) INTEGER
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	75 * IHIZ (input) INTEGER
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	76 * Specify the rows of Z to which transformations must be
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	77 * applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	78 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	79 * Z (input/output) COMPLEX array of size (LDZ,IHI)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	80 * If WANTZ = .TRUE., then the QR Sweep unitary
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	81 * similarity transformation is accumulated into
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	82 * Z(ILOZ:IHIZ,ILO:IHI) from the right.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	83 * If WANTZ = .FALSE., then Z is unreferenced.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	84 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	85 * LDZ (input) integer scalar
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	86 * LDA is the leading dimension of Z just as declared in
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	87 * the calling procedure. LDZ.GE.N.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	88 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	89 * V (workspace) COMPLEX array of size (LDV,NSHFTS/2)
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 * LDV (input) integer scalar
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	92 * LDV is the leading dimension of V as declared in the
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	93 * calling procedure. LDV.GE.3.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	94 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	95 * U (workspace) COMPLEX array of size
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	96 * (LDU,3*NSHFTS-3)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	97 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	98 * LDU (input) integer scalar
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	99 * LDU is the leading dimension of U just as declared in the
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	100 * in the calling subroutine. LDU.GE.3*NSHFTS-3.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	101 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	102 * NH (input) integer scalar
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	103 * NH is the number of columns in array WH available for
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	104 * workspace. NH.GE.1.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	105 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	106 * WH (workspace) COMPLEX array of size (LDWH,NH)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	107 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	108 * LDWH (input) integer scalar
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	109 * Leading dimension of WH just as declared in the
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	110 * calling procedure. LDWH.GE.3*NSHFTS-3.
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 * NV (input) integer scalar
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	113 * NV is the number of rows in WV agailable for workspace.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	114 * NV.GE.1.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	115 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	116 * WV (workspace) COMPLEX array of size
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	117 * (LDWV,3*NSHFTS-3)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	118 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	119 * LDWV (input) integer scalar
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	120 * LDWV is the leading dimension of WV as declared in the
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	121 * in the calling subroutine. LDWV.GE.NV.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	122 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	123 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	124 * ================================================================
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	125 * Based on contributions by
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	126 * Karen Braman and Ralph Byers, Department of Mathematics,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	127 * University of Kansas, USA
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	128 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	129 * ============================================================
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	130 * Reference:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	131 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	132 * 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	133 * Algorithm Part I: Maintaining Well Focused Shifts, and
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	134 * Level 3 Performance, SIAM Journal of Matrix Analysis,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	135 * volume 23, pages 929--947, 2002.
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 * .. Parameters ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	139 COMPLEX ZERO, ONE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	140 PARAMETER ( ZERO = ( 0.0e0, 0.0e0 ),
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	141 $ ONE = ( 1.0e0, 0.0e0 ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	142 REAL RZERO, RONE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	143 PARAMETER ( RZERO = 0.0e0, RONE = 1.0e0 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	144 * ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	145 * .. Local Scalars ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	146 COMPLEX ALPHA, BETA, CDUM, REFSUM
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	147 REAL H11, H12, H21, H22, SAFMAX, SAFMIN, SCL,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	148 $ SMLNUM, TST1, TST2, ULP
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	149 INTEGER I2, I4, INCOL, J, J2, J4, JBOT, JCOL, JLEN,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	150 $ JROW, JTOP, K, K1, KDU, KMS, KNZ, KRCOL, KZS,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	151 $ M, M22, MBOT, MEND, MSTART, MTOP, NBMPS, NDCOL,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	152 $ NS, NU
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	153 LOGICAL ACCUM, BLK22, BMP22
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	154 * ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	155 * .. External Functions ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	156 REAL SLAMCH
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	157 EXTERNAL SLAMCH
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	158 * ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	159 * .. Intrinsic Functions ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	160 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	161 INTRINSIC ABS, AIMAG, CONJG, MAX, MIN, MOD, REAL
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	162 * ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	163 * .. Local Arrays ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	164 COMPLEX VT( 3 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	165 * ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	166 * .. External Subroutines ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	167 EXTERNAL CGEMM, CLACPY, CLAQR1, CLARFG, CLASET, CTRMM,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	168 $ SLABAD
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	169 * ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	170 * .. Statement Functions ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	171 REAL CABS1
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	172 * ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	173 * .. Statement Function definitions ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	174 CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	175 * ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	176 * .. Executable Statements ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	177 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	178 * ==== If there are no shifts, then there is nothing to do. ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	179 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	180 IF( NSHFTS.LT.2 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	181 $ RETURN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	182 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	183 * ==== If the active block is empty or 1-by-1, then there
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	184 * . is nothing to do. ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	185 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	186 IF( KTOP.GE.KBOT )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	187 $ RETURN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	188 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	189 * ==== NSHFTS is supposed to be even, but if is odd,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	190 * . then simply reduce it by one. ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	191 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	192 NS = NSHFTS - MOD( NSHFTS, 2 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	193 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	194 * ==== Machine constants for deflation ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	195 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	196 SAFMIN = SLAMCH( 'SAFE MINIMUM' )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	197 SAFMAX = RONE / SAFMIN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	198 CALL SLABAD( SAFMIN, SAFMAX )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	199 ULP = SLAMCH( 'PRECISION' )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	200 SMLNUM = SAFMIN*( REAL( N ) / ULP )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	201 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	202 * ==== Use accumulated reflections to update far-from-diagonal
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	203 * . entries ? ====
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 ACCUM = ( KACC22.EQ.1 ) .OR. ( KACC22.EQ.2 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	206 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	207 * ==== If so, exploit the 2-by-2 block structure? ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	208 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	209 BLK22 = ( NS.GT.2 ) .AND. ( KACC22.EQ.2 )
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 * ==== clear trash ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	212 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	213 IF( KTOP+2.LE.KBOT )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	214 $ H( KTOP+2, KTOP ) = ZERO
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 * ==== NBMPS = number of 2-shift bulges in the chain ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	217 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	218 NBMPS = NS / 2
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	219 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	220 * ==== KDU = width of slab ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	221 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	222 KDU = 6*NBMPS - 3
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	223 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	224 * ==== Create and chase chains of NBMPS bulges ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	225 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	226 DO 210 INCOL = 3( 1-NBMPS ) + KTOP - 1, KBOT - 2, 3NBMPS - 2
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	227 NDCOL = INCOL + KDU
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	228 IF( ACCUM )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	229 $ CALL CLASET( 'ALL', KDU, KDU, ZERO, ONE, U, LDU )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	230 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	231 * ==== Near-the-diagonal bulge chase. The following loop
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	232 * . performs the near-the-diagonal part of a small bulge
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	233 * . multi-shift QR sweep. Each 6*NBMPS-2 column diagonal
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	234 * . chunk extends from column INCOL to column NDCOL
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	235 * . (including both column INCOL and column NDCOL). The
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	236 * . following loop chases a 3*NBMPS column long chain of
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	237 * . NBMPS bulges 3*NBMPS-2 columns to the right. (INCOL
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	238 * . may be less than KTOP and and NDCOL may be greater than
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	239 * . KBOT indicating phantom columns from which to chase
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	240 * . bulges before they are actually introduced or to which
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	241 * . to chase bulges beyond column KBOT.) ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	242 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	243 DO 140 KRCOL = INCOL, MIN( INCOL+3*NBMPS-3, KBOT-2 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	244 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	245 * ==== Bulges number MTOP to MBOT are active double implicit
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	246 * . shift bulges. There may or may not also be small
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	247 * . 2-by-2 bulge, if there is room. The inactive bulges
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	248 * . (if any) must wait until the active bulges have moved
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	249 * . down the diagonal to make room. The phantom matrix
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	250 * . paradigm described above helps keep track. ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	251 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	252 MTOP = MAX( 1, ( ( KTOP-1 )-KRCOL+2 ) / 3+1 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	253 MBOT = MIN( NBMPS, ( KBOT-KRCOL ) / 3 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	254 M22 = MBOT + 1
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	255 BMP22 = ( MBOT.LT.NBMPS ) .AND. ( KRCOL+3*( M22-1 ) ).EQ.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	256 $ ( KBOT-2 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	257 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	258 * ==== Generate reflections to chase the chain right
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	259 * . one column. (The minimum value of K is KTOP-1.) ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	260 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	261 DO 10 M = MTOP, MBOT
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	262 K = KRCOL + 3*( M-1 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	263 IF( K.EQ.KTOP-1 ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	264 CALL CLAQR1( 3, H( KTOP, KTOP ), LDH, S( 2*M-1 ),
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	265 $ S( 2*M ), V( 1, M ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	266 ALPHA = V( 1, M )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	267 CALL CLARFG( 3, ALPHA, V( 2, M ), 1, V( 1, M ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	268 ELSE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	269 BETA = H( K+1, K )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	270 V( 2, M ) = H( K+2, K )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	271 V( 3, M ) = H( K+3, K )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	272 CALL CLARFG( 3, BETA, V( 2, M ), 1, V( 1, M ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	273 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	274 * ==== A Bulge may collapse because of vigilant
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	275 * . deflation or destructive underflow. (The
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	276 * . initial bulge is always collapsed.) Use
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	277 * . the two-small-subdiagonals trick to try
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	278 * . to get it started again. If V(2,M).NE.0 and
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	279 * . V(3,M) = H(K+3,K+1) = H(K+3,K+2) = 0, then
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	280 * . this bulge is collapsing into a zero
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	281 * . subdiagonal. It will be restarted next
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	282 * . trip through the loop.)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	283 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	284 IF( V( 1, M ).NE.ZERO .AND.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	285 $ ( V( 3, M ).NE.ZERO .OR. ( H( K+3,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	286 $ K+1 ).EQ.ZERO .AND. H( K+3, K+2 ).EQ.ZERO ) ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	287 $ THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	288 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	289 * ==== Typical case: not collapsed (yet). ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	290 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	291 H( K+1, K ) = BETA
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	292 H( K+2, K ) = ZERO
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	293 H( K+3, K ) = ZERO
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	294 ELSE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	295 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	296 * ==== Atypical case: collapsed. Attempt to
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	297 * . reintroduce ignoring H(K+1,K). If the
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	298 * . fill resulting from the new reflector
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	299 * . is too large, then abandon it.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	300 * . Otherwise, use the new one. ====
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 CALL CLAQR1( 3, H( K+1, K+1 ), LDH, S( 2*M-1 ),
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	303 $ S( 2*M ), VT )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	304 SCL = CABS1( VT( 1 ) ) + CABS1( VT( 2 ) ) +
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	305 $ CABS1( VT( 3 ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	306 IF( SCL.NE.RZERO ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	307 VT( 1 ) = VT( 1 ) / SCL
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	308 VT( 2 ) = VT( 2 ) / SCL
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	309 VT( 3 ) = VT( 3 ) / SCL
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	310 END IF
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 * ==== The following is the traditional and
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	313 * . conservative two-small-subdiagonals
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	314 * . test. ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	315 * .
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	316 IF( CABS1( H( K+1, K ) )*
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	317 $ ( CABS1( VT( 2 ) )+CABS1( VT( 3 ) ) ).GT.ULP*
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	318 $ CABS1( VT( 1 ) )*( CABS1( H( K,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	319 $ K ) )+CABS1( H( K+1, K+1 ) )+CABS1( H( K+2,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	320 $ K+2 ) ) ) ) THEN
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 * ==== Starting a new bulge here would
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	323 * . create non-negligible fill. If
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	324 * . the old reflector is diagonal (only
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	325 * . possible with underflows), then
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	326 * . change it to I. Otherwise, use
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	327 * . it with trepidation. ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	328 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	329 IF( V( 2, M ).EQ.ZERO .AND. V( 3, M ).EQ.ZERO )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	330 $ THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	331 V( 1, M ) = ZERO
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	332 ELSE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	333 H( K+1, K ) = BETA
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	334 H( K+2, K ) = ZERO
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	335 H( K+3, K ) = ZERO
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	336 END IF
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	337 ELSE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	338 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	339 * ==== Stating a new bulge here would
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	340 * . create only negligible fill.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	341 * . Replace the old reflector with
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	342 * . the new one. ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	343 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	344 ALPHA = VT( 1 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	345 CALL CLARFG( 3, ALPHA, VT( 2 ), 1, VT( 1 ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	346 REFSUM = H( K+1, K ) +
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	347 $ H( K+2, K )*CONJG( VT( 2 ) ) +
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	348 $ H( K+3, K )*CONJG( VT( 3 ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	349 H( K+1, K ) = H( K+1, K ) -
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	350 $ CONJG( VT( 1 ) )*REFSUM
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	351 H( K+2, K ) = ZERO
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	352 H( K+3, K ) = ZERO
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	353 V( 1, M ) = VT( 1 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	354 V( 2, M ) = VT( 2 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	355 V( 3, M ) = VT( 3 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	356 END IF
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	357 END IF
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	358 END IF
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	359 10 CONTINUE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	360 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	361 * ==== Generate a 2-by-2 reflection, if needed. ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	362 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	363 K = KRCOL + 3*( M22-1 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	364 IF( BMP22 ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	365 IF( K.EQ.KTOP-1 ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	366 CALL CLAQR1( 2, H( K+1, K+1 ), LDH, S( 2*M22-1 ),
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	367 $ S( 2*M22 ), V( 1, M22 ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	368 BETA = V( 1, M22 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	369 CALL CLARFG( 2, BETA, V( 2, M22 ), 1, V( 1, M22 ) )
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 BETA = H( K+1, K )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	372 V( 2, M22 ) = H( K+2, K )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	373 CALL CLARFG( 2, BETA, V( 2, M22 ), 1, V( 1, M22 ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	374 H( K+1, K ) = BETA
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	375 H( K+2, K ) = ZERO
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 * ==== Initialize V(1,M22) here to avoid possible undefined
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	380 * . variable problems later. ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	381 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	382 V( 1, M22 ) = ZERO
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	383 END IF
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	384 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	385 * ==== Multiply H by reflections from the left ====
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( ACCUM ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	388 JBOT = MIN( NDCOL, KBOT )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	389 ELSE IF( WANTT ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	390 JBOT = N
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	391 ELSE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	392 JBOT = KBOT
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 DO 30 J = MAX( KTOP, KRCOL ), JBOT
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	395 MEND = MIN( MBOT, ( J-KRCOL+2 ) / 3 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	396 DO 20 M = MTOP, MEND
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	397 K = KRCOL + 3*( M-1 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	398 REFSUM = CONJG( V( 1, M ) )*
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	399 $ ( H( K+1, J )+CONJG( V( 2, M ) )*H( K+2, J )+
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	400 $ CONJG( V( 3, M ) )*H( K+3, J ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	401 H( K+1, J ) = H( K+1, J ) - REFSUM
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	402 H( K+2, J ) = H( K+2, J ) - REFSUM*V( 2, M )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	403 H( K+3, J ) = H( K+3, J ) - REFSUM*V( 3, M )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	404 20 CONTINUE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	405 30 CONTINUE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	406 IF( BMP22 ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	407 K = KRCOL + 3*( M22-1 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	408 DO 40 J = MAX( K+1, KTOP ), JBOT
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	409 REFSUM = CONJG( V( 1, M22 ) )*
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	410 $ ( H( K+1, J )+CONJG( V( 2, M22 ) )*
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	411 $ H( K+2, J ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	412 H( K+1, J ) = H( K+1, J ) - REFSUM
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	413 H( K+2, J ) = H( K+2, J ) - REFSUM*V( 2, M22 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	414 40 CONTINUE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	415 END IF
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	416 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	417 * ==== Multiply H by reflections from the right.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	418 * . Delay filling in the last row until the
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	419 * . vigilant deflation check is complete. ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	420 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	421 IF( ACCUM ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	422 JTOP = MAX( KTOP, INCOL )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	423 ELSE IF( WANTT ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	424 JTOP = 1
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	425 ELSE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	426 JTOP = KTOP
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	427 END IF
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	428 DO 80 M = MTOP, MBOT
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	429 IF( V( 1, M ).NE.ZERO ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	430 K = KRCOL + 3*( M-1 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	431 DO 50 J = JTOP, MIN( KBOT, K+3 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	432 REFSUM = V( 1, M )( H( J, K+1 )+V( 2, M )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	433 $ H( J, K+2 )+V( 3, M )*H( J, K+3 ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	434 H( J, K+1 ) = H( J, K+1 ) - REFSUM
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	435 H( J, K+2 ) = H( J, K+2 ) -
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	436 $ REFSUM*CONJG( V( 2, M ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	437 H( J, K+3 ) = H( J, K+3 ) -
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	438 $ REFSUM*CONJG( V( 3, M ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	439 50 CONTINUE
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 IF( ACCUM ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	442 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	443 * ==== Accumulate U. (If necessary, update Z later
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	444 * . with with an efficient matrix-matrix
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	445 * . multiply.) ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	446 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	447 KMS = K - INCOL
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	448 DO 60 J = MAX( 1, KTOP-INCOL ), KDU
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	449 REFSUM = V( 1, M )( U( J, KMS+1 )+V( 2, M )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	450 $ U( J, KMS+2 )+V( 3, M )*U( J, KMS+3 ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	451 U( J, KMS+1 ) = U( J, KMS+1 ) - REFSUM
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	452 U( J, KMS+2 ) = U( J, KMS+2 ) -
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	453 $ REFSUM*CONJG( V( 2, M ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	454 U( J, KMS+3 ) = U( J, KMS+3 ) -
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	455 $ REFSUM*CONJG( V( 3, M ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	456 60 CONTINUE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	457 ELSE IF( WANTZ ) THEN
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 * ==== U is not accumulated, so update Z
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	460 * . now by multiplying by reflections
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	461 * . from the right. ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	462 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	463 DO 70 J = ILOZ, IHIZ
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	464 REFSUM = V( 1, M )( Z( J, K+1 )+V( 2, M )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	465 $ Z( J, K+2 )+V( 3, M )*Z( J, K+3 ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	466 Z( J, K+1 ) = Z( J, K+1 ) - REFSUM
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	467 Z( J, K+2 ) = Z( J, K+2 ) -
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	468 $ REFSUM*CONJG( V( 2, M ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	469 Z( J, K+3 ) = Z( J, K+3 ) -
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	470 $ REFSUM*CONJG( V( 3, M ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	471 70 CONTINUE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	472 END IF
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	473 END IF
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	474 80 CONTINUE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	475 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	476 * ==== Special case: 2-by-2 reflection (if needed) ====
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	477 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	478 K = KRCOL + 3*( M22-1 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	479 IF( BMP22 .AND. ( V( 1, M22 ).NE.ZERO ) ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	480 DO 90 J = JTOP, MIN( KBOT, K+3 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	481 REFSUM = V( 1, M22 )( H( J, K+1 )+V( 2, M22 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	482 $ H( J, K+2 ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	483 H( J, K+1 ) = H( J, K+1 ) - REFSUM
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	484 H( J, K+2 ) = H( J, K+2 ) -
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	485 $ REFSUM*CONJG( V( 2, M22 ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	486 90 CONTINUE
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( ACCUM ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	489 KMS = K - INCOL
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	490 DO 100 J = MAX( 1, KTOP-INCOL ), KDU
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	491 REFSUM = V( 1, M22 )( U( J, KMS+1 )+V( 2, M22 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	492 $ U( J, KMS+2 ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	493 U( J, KMS+1 ) = U( J, KMS+1 ) - REFSUM
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	494 U( J, KMS+2 ) = U( J, KMS+2 ) -
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	495 $ REFSUM*CONJG( V( 2, M22 ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	496 100 CONTINUE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	497 ELSE IF( WANTZ ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	498 DO 110 J = ILOZ, IHIZ
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	499 REFSUM = V( 1, M22 )( Z( J, K+1 )+V( 2, M22 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	500 $ Z( J, K+2 ) )

7789

82be108cc558 First attempt at single precision tyeps