octave-nkf: libcruft/lapack/cbdsqr.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 CBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	2 $ LDU, C, LDC, RWORK, 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 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 CHARACTER UPLO
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	10 INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU
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 REAL D( * ), E( * ), RWORK( * )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	14 COMPLEX C( LDC, * ), U( LDU, * ), VT( LDVT, * )
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 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	17 * Purpose
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 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	20 * CBDSQR computes the singular values and, optionally, the right and/or
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	21 * left singular vectors from the singular value decomposition (SVD) of
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	22 * a real N-by-N (upper or lower) bidiagonal matrix B using the implicit
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	23 * zero-shift QR algorithm. The SVD of B has the form
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	24 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	25 * B = Q * S * P**H
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	26 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	27 * where S is the diagonal matrix of singular values, Q is an orthogonal
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	28 * matrix of left singular vectors, and P is an orthogonal matrix of
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	29 * right singular vectors. If left singular vectors are requested, this
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	30 * subroutine actually returns U*Q instead of Q, and, if right singular
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	31 * vectors are requested, this subroutine returns P*HVT instead of
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	32 * P**H, for given complex input matrices U and VT. When U and VT are
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	33 * the unitary matrices that reduce a general matrix A to bidiagonal
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	34 * form: A = UBVT, as computed by CGEBRD, then
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 * A = (UQ) S * (P*HVT)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	37 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	38 * is the SVD of A. Optionally, the subroutine may also compute Q*HC
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	39 * for a given complex input matrix C.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	40 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	41 * See "Computing Small Singular Values of Bidiagonal Matrices With
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	42 * Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	43 * LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	44 * no. 5, pp. 873-912, Sept 1990) and
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	45 * "Accurate singular values and differential qd algorithms," by
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	46 * B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	47 * Department, University of California at Berkeley, July 1992
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	48 * for a detailed description of the algorithm.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	49 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	50 * Arguments
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	51 * =========
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	52 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	53 * UPLO (input) CHARACTER*1
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	54 * = 'U': B is upper bidiagonal;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	55 * = 'L': B is lower bidiagonal.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	56 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	57 * N (input) INTEGER
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	58 * The order of the matrix B. N >= 0.
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 * NCVT (input) INTEGER
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	61 * The number of columns of the matrix VT. NCVT >= 0.
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 * NRU (input) INTEGER
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	64 * The number of rows of the matrix U. NRU >= 0.
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 * NCC (input) INTEGER
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	67 * The number of columns of the matrix C. NCC >= 0.
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 * D (input/output) REAL array, dimension (N)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	70 * On entry, the n diagonal elements of the bidiagonal matrix B.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	71 * On exit, if INFO=0, the singular values of B in decreasing
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	72 * order.
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 * E (input/output) REAL array, dimension (N-1)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	75 * On entry, the N-1 offdiagonal elements of the bidiagonal
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	76 * matrix B.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	77 * On exit, if INFO = 0, E is destroyed; if INFO > 0, D and E
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	78 * will contain the diagonal and superdiagonal elements of a
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	79 * bidiagonal matrix orthogonally equivalent to the one given
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	80 * as input.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	81 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	82 * VT (input/output) COMPLEX array, dimension (LDVT, NCVT)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	83 * On entry, an N-by-NCVT matrix VT.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	84 * On exit, VT is overwritten by P*H VT.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	85 * Not referenced if NCVT = 0.
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 * LDVT (input) INTEGER
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	88 * The leading dimension of the array VT.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	89 * LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0.
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 * U (input/output) COMPLEX array, dimension (LDU, N)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	92 * On entry, an NRU-by-N matrix U.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	93 * On exit, U is overwritten by U * Q.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	94 * Not referenced if NRU = 0.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	95 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	96 * LDU (input) INTEGER
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	97 * The leading dimension of the array U. LDU >= max(1,NRU).
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	98 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	99 * C (input/output) COMPLEX array, dimension (LDC, NCC)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	100 * On entry, an N-by-NCC matrix C.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	101 * On exit, C is overwritten by Q*H C.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	102 * Not referenced if NCC = 0.
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 * LDC (input) INTEGER
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	105 * The leading dimension of the array C.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	106 * LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0.
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 * RWORK (workspace) REAL array, dimension (2*N)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	109 * if NCVT = NRU = NCC = 0, (max(1, 4*N-4)) otherwise
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	110 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	111 * INFO (output) INTEGER
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	112 * = 0: successful exit
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	113 * < 0: If INFO = -i, the i-th argument had an illegal value
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	114 * > 0: the algorithm did not converge; D and E contain the
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	115 * elements of a bidiagonal matrix which is orthogonally
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	116 * similar to the input matrix B; if INFO = i, i
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	117 * elements of E have not converged to zero.
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 * Internal Parameters
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	120 * ===================
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 * TOLMUL REAL, default = max(10,min(100,EPS**(-1/8)))
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	123 * TOLMUL controls the convergence criterion of the QR loop.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	124 * If it is positive, TOLMUL*EPS is the desired relative
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	125 * precision in the computed singular values.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	126 * If it is negative, abs(TOLMULEPSsigma_max) is the
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	127 * desired absolute accuracy in the computed singular
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	128 * values (corresponds to relative accuracy
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	129 * abs(TOLMUL*EPS) in the largest singular value.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	130 * abs(TOLMUL) should be between 1 and 1/EPS, and preferably
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	131 * between 10 (for fast convergence) and .1/EPS
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	132 * (for there to be some accuracy in the results).
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	133 * Default is to lose at either one eighth or 2 of the
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	134 * available decimal digits in each computed singular value
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	135 * (whichever is smaller).
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 * MAXITR INTEGER, default = 6
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	138 * MAXITR controls the maximum number of passes of the
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	139 * algorithm through its inner loop. The algorithms stops
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	140 * (and so fails to converge) if the number of passes
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	141 * through the inner loop exceeds MAXITRN*2.
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 * =====================================================================
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 * .. Parameters ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	146 REAL ZERO
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	147 PARAMETER ( ZERO = 0.0E0 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	148 REAL ONE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	149 PARAMETER ( ONE = 1.0E0 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	150 REAL NEGONE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	151 PARAMETER ( NEGONE = -1.0E0 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	152 REAL HNDRTH
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	153 PARAMETER ( HNDRTH = 0.01E0 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	154 REAL TEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	155 PARAMETER ( TEN = 10.0E0 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	156 REAL HNDRD
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	157 PARAMETER ( HNDRD = 100.0E0 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	158 REAL MEIGTH
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	159 PARAMETER ( MEIGTH = -0.125E0 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	160 INTEGER MAXITR
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	161 PARAMETER ( MAXITR = 6 )
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 Scalars ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	164 LOGICAL LOWER, ROTATE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	165 INTEGER I, IDIR, ISUB, ITER, J, LL, LLL, M, MAXIT, NM1,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	166 $ NM12, NM13, OLDLL, OLDM
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	167 REAL ABSE, ABSS, COSL, COSR, CS, EPS, F, G, H, MU,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	168 $ OLDCS, OLDSN, R, SHIFT, SIGMN, SIGMX, SINL,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	169 $ SINR, SLL, SMAX, SMIN, SMINL, SMINOA,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	170 $ SN, THRESH, TOL, TOLMUL, UNFL
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	171 * ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	172 * .. External Functions ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	173 LOGICAL LSAME
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	174 REAL SLAMCH
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	175 EXTERNAL LSAME, SLAMCH
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	176 * ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	177 * .. External Subroutines ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	178 EXTERNAL CLASR, CSROT, CSSCAL, CSWAP, SLARTG, SLAS2,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	179 $ SLASQ1, SLASV2, XERBLA
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 * .. Intrinsic Functions ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	182 INTRINSIC ABS, MAX, MIN, REAL, SIGN, SQRT
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	183 * ..
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	184 * .. Executable Statements ..
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 * Test the input parameters.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	187 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	188 INFO = 0
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	189 LOWER = LSAME( UPLO, 'L' )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	190 IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LOWER ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	191 INFO = -1
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	192 ELSE IF( N.LT.0 ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	193 INFO = -2
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	194 ELSE IF( NCVT.LT.0 ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	195 INFO = -3
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	196 ELSE IF( NRU.LT.0 ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	197 INFO = -4
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	198 ELSE IF( NCC.LT.0 ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	199 INFO = -5
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	200 ELSE IF( ( NCVT.EQ.0 .AND. LDVT.LT.1 ) .OR.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	201 $ ( NCVT.GT.0 .AND. LDVT.LT.MAX( 1, N ) ) ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	202 INFO = -9
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	203 ELSE IF( LDU.LT.MAX( 1, NRU ) ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	204 INFO = -11
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	205 ELSE IF( ( NCC.EQ.0 .AND. LDC.LT.1 ) .OR.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	206 $ ( NCC.GT.0 .AND. LDC.LT.MAX( 1, N ) ) ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	207 INFO = -13
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	208 END IF
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	209 IF( INFO.NE.0 ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	210 CALL XERBLA( 'CBDSQR', -INFO )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	211 RETURN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	212 END IF
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	213 IF( N.EQ.0 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	214 $ RETURN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	215 IF( N.EQ.1 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	216 $ GO TO 160
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 * ROTATE is true if any singular vectors desired, false otherwise
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 ROTATE = ( NCVT.GT.0 ) .OR. ( NRU.GT.0 ) .OR. ( NCC.GT.0 )
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 * If no singular vectors desired, use qd algorithm
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 IF( .NOT.ROTATE ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	225 CALL SLASQ1( N, D, E, RWORK, INFO )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	226 RETURN
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 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	229 NM1 = N - 1
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	230 NM12 = NM1 + NM1
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	231 NM13 = NM12 + NM1
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	232 IDIR = 0
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 * Get machine constants
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 EPS = SLAMCH( 'Epsilon' )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	237 UNFL = SLAMCH( 'Safe minimum' )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	238 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	239 * If matrix lower bidiagonal, rotate to be upper bidiagonal
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	240 * by applying Givens rotations on the left
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	241 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	242 IF( LOWER ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	243 DO 10 I = 1, N - 1
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	244 CALL SLARTG( D( I ), E( I ), CS, SN, R )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	245 D( I ) = R
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	246 E( I ) = SN*D( I+1 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	247 D( I+1 ) = CS*D( I+1 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	248 RWORK( I ) = CS
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	249 RWORK( NM1+I ) = SN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	250 10 CONTINUE
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 * Update singular vectors if desired
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	253 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	254 IF( NRU.GT.0 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	255 $ CALL CLASR( 'R', 'V', 'F', NRU, N, RWORK( 1 ), RWORK( N ),
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	256 $ U, LDU )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	257 IF( NCC.GT.0 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	258 $ CALL CLASR( 'L', 'V', 'F', N, NCC, RWORK( 1 ), RWORK( N ),
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	259 $ C, LDC )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	260 END IF
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	261 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	262 * Compute singular values to relative accuracy TOL
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	263 * (By setting TOL to be negative, algorithm will compute
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	264 * singular values to absolute accuracy ABS(TOL)*norm(input matrix))
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	265 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	266 TOLMUL = MAX( TEN, MIN( HNDRD, EPS**MEIGTH ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	267 TOL = TOLMUL*EPS
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 * Compute approximate maximum, minimum singular values
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	270 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	271 SMAX = ZERO
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	272 DO 20 I = 1, N
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	273 SMAX = MAX( SMAX, ABS( D( I ) ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	274 20 CONTINUE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	275 DO 30 I = 1, N - 1
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	276 SMAX = MAX( SMAX, ABS( E( I ) ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	277 30 CONTINUE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	278 SMINL = ZERO
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	279 IF( TOL.GE.ZERO ) THEN
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 * Relative accuracy desired
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 SMINOA = ABS( D( 1 ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	284 IF( SMINOA.EQ.ZERO )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	285 $ GO TO 50
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	286 MU = SMINOA
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	287 DO 40 I = 2, N
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	288 MU = ABS( D( I ) )*( MU / ( MU+ABS( E( I-1 ) ) ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	289 SMINOA = MIN( SMINOA, MU )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	290 IF( SMINOA.EQ.ZERO )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	291 $ GO TO 50
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	292 40 CONTINUE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	293 50 CONTINUE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	294 SMINOA = SMINOA / SQRT( REAL( N ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	295 THRESH = MAX( TOLSMINOA, MAXITRNNUNFL )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	296 ELSE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	297 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	298 * Absolute accuracy desired
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	299 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	300 THRESH = MAX( ABS( TOL )SMAX, MAXITRNNUNFL )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	301 END IF
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	302 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	303 * Prepare for main iteration loop for the singular values
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	304 * (MAXIT is the maximum number of passes through the inner
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	305 * loop permitted before nonconvergence signalled.)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	306 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	307 MAXIT = MAXITRNN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	308 ITER = 0
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	309 OLDLL = -1
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	310 OLDM = -1
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 * M points to last element of unconverged part of matrix
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 M = N
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 * Begin main iteration loop
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	317 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	318 60 CONTINUE
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 * Check for convergence or exceeding iteration count
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 IF( M.LE.1 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	323 $ GO TO 160
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	324 IF( ITER.GT.MAXIT )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	325 $ GO TO 200
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	326 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	327 * Find diagonal block of matrix to work on
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( TOL.LT.ZERO .AND. ABS( D( M ) ).LE.THRESH )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	330 $ D( M ) = ZERO
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	331 SMAX = ABS( D( M ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	332 SMIN = SMAX
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	333 DO 70 LLL = 1, M - 1
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	334 LL = M - LLL
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	335 ABSS = ABS( D( LL ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	336 ABSE = ABS( E( LL ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	337 IF( TOL.LT.ZERO .AND. ABSS.LE.THRESH )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	338 $ D( LL ) = ZERO
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	339 IF( ABSE.LE.THRESH )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	340 $ GO TO 80
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	341 SMIN = MIN( SMIN, ABSS )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	342 SMAX = MAX( SMAX, ABSS, ABSE )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	343 70 CONTINUE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	344 LL = 0
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	345 GO TO 90
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	346 80 CONTINUE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	347 E( LL ) = ZERO
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	348 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	349 * Matrix splits since E(LL) = 0
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 IF( LL.EQ.M-1 ) THEN
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 * Convergence of bottom singular value, return to top of loop
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	354 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	355 M = M - 1
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	356 GO TO 60
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 90 CONTINUE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	359 LL = LL + 1
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 * E(LL) through E(M-1) are nonzero, E(LL-1) is zero
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 IF( LL.EQ.M-1 ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	364 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	365 * 2 by 2 block, handle separately
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	366 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	367 CALL SLASV2( D( M-1 ), E( M-1 ), D( M ), SIGMN, SIGMX, SINR,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	368 $ COSR, SINL, COSL )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	369 D( M-1 ) = SIGMX
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	370 E( M-1 ) = ZERO
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	371 D( M ) = SIGMN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	372 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	373 * Compute singular vectors, if desired
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	374 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	375 IF( NCVT.GT.0 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	376 $ CALL CSROT( NCVT, VT( M-1, 1 ), LDVT, VT( M, 1 ), LDVT,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	377 $ COSR, SINR )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	378 IF( NRU.GT.0 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	379 $ CALL CSROT( NRU, U( 1, M-1 ), 1, U( 1, M ), 1, COSL, SINL )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	380 IF( NCC.GT.0 )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	381 $ CALL CSROT( NCC, C( M-1, 1 ), LDC, C( M, 1 ), LDC, COSL,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	382 $ SINL )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	383 M = M - 2
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	384 GO TO 60
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	385 END IF
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 working on new submatrix, choose shift direction
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	388 * (from larger end diagonal element towards smaller)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	389 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	390 IF( LL.GT.OLDM .OR. M.LT.OLDLL ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	391 IF( ABS( D( LL ) ).GE.ABS( D( M ) ) ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	392 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	393 * Chase bulge from top (big end) to bottom (small end)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	394 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	395 IDIR = 1
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	396 ELSE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	397 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	398 * Chase bulge from bottom (big end) to top (small end)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	399 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	400 IDIR = 2
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	401 END IF
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	402 END IF
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	403 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	404 * Apply convergence tests
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 IF( IDIR.EQ.1 ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	407 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	408 * Run convergence test in forward direction
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	409 * First apply standard test to bottom of matrix
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	410 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	411 IF( ABS( E( M-1 ) ).LE.ABS( TOL )*ABS( D( M ) ) .OR.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	412 $ ( TOL.LT.ZERO .AND. ABS( E( M-1 ) ).LE.THRESH ) ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	413 E( M-1 ) = ZERO
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	414 GO TO 60
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 IF( TOL.GE.ZERO ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	418 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	419 * If relative accuracy desired,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	420 * apply convergence criterion forward
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 MU = ABS( D( LL ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	423 SMINL = MU
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	424 DO 100 LLL = LL, M - 1
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	425 IF( ABS( E( LLL ) ).LE.TOL*MU ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	426 E( LLL ) = ZERO
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	427 GO TO 60
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	428 END IF
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	429 MU = ABS( D( LLL+1 ) )*( MU / ( MU+ABS( E( LLL ) ) ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	430 SMINL = MIN( SMINL, MU )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	431 100 CONTINUE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	432 END IF
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 ELSE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	435 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	436 * Run convergence test in backward direction
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	437 * First apply standard test to top of matrix
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	438 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	439 IF( ABS( E( LL ) ).LE.ABS( TOL )*ABS( D( LL ) ) .OR.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	440 $ ( TOL.LT.ZERO .AND. ABS( E( LL ) ).LE.THRESH ) ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	441 E( LL ) = ZERO
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	442 GO TO 60
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	443 END IF
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	444 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	445 IF( TOL.GE.ZERO ) THEN
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 * If relative accuracy desired,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	448 * apply convergence criterion backward
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	449 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	450 MU = ABS( D( M ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	451 SMINL = MU
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	452 DO 110 LLL = M - 1, LL, -1
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	453 IF( ABS( E( LLL ) ).LE.TOL*MU ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	454 E( LLL ) = ZERO
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	455 GO TO 60
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	456 END IF
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	457 MU = ABS( D( LLL ) )*( MU / ( MU+ABS( E( LLL ) ) ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	458 SMINL = MIN( SMINL, MU )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	459 110 CONTINUE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	460 END IF
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	461 END IF
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	462 OLDLL = LL
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	463 OLDM = M
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 * Compute shift. First, test if shifting would ruin relative
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	466 * accuracy, and if so set the shift to zero.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	467 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	468 IF( TOL.GE.ZERO .AND. NTOL( SMINL / SMAX ).LE.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	469 $ MAX( EPS, HNDRTH*TOL ) ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	470 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	471 * Use a zero shift to avoid loss of relative accuracy
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	472 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	473 SHIFT = ZERO
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 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	476 * Compute the shift from 2-by-2 block at end of matrix
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 IF( IDIR.EQ.1 ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	479 SLL = ABS( D( LL ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	480 CALL SLAS2( D( M-1 ), E( M-1 ), D( M ), SHIFT, R )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	481 ELSE
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	482 SLL = ABS( D( M ) )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	483 CALL SLAS2( D( LL ), E( LL ), D( LL+1 ), SHIFT, R )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	484 END IF
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	485 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	486 * Test if shift negligible, and if so set to zero
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( SLL.GT.ZERO ) THEN
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	489 IF( ( SHIFT / SLL )**2.LT.EPS )
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	490 $ SHIFT = ZERO
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	491 END IF
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	492 END IF
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	493 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	494 * Increment iteration count
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	495 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	496 ITER = ITER + M - LL
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	497 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	498 * If SHIFT = 0, do simplified QR iteration
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	499 *
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: diff changeset	500 IF( SHIFT.EQ.ZERO ) THEN

7789

82be108cc558 First attempt at single precision tyeps