Mercurial > octave-nkf

annotate libcruft/lapack/dlasd7.f @ 7072:b48d486f641d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2007-10-26 15:52:57 by jwe]
author jwe
date Fri, 26 Oct 2007 15:52:58 +0000
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
7072
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
1 SUBROUTINE DLASD7( ICOMPQ, NL, NR, SQRE, K, D, Z, ZW, VF, VFW, VL,
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
2 $ VLW, ALPHA, BETA, DSIGMA, IDX, IDXP, IDXQ,
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
3 $ PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM,
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
4 $ C, S, INFO )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
5 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
6 * -- LAPACK auxiliary routine (version 3.1) --
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
7 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
8 * November 2006
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
9 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
10 * .. Scalar Arguments ..
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
11 INTEGER GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL,
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
12 $ NR, SQRE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
13 DOUBLE PRECISION ALPHA, BETA, C, S
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
14 * ..
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
15 * .. Array Arguments ..
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
16 INTEGER GIVCOL( LDGCOL, * ), IDX( * ), IDXP( * ),
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
17 $ IDXQ( * ), PERM( * )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
18 DOUBLE PRECISION D( * ), DSIGMA( * ), GIVNUM( LDGNUM, * ),
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
19 $ VF( * ), VFW( * ), VL( * ), VLW( * ), Z( * ),
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
20 $ ZW( * )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
21 * ..
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
22 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
23 * Purpose
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
24 * =======
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
25 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
26 * DLASD7 merges the two sets of singular values together into a single
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
27 * sorted set. Then it tries to deflate the size of the problem. There
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
28 * are two ways in which deflation can occur: when two or more singular
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
29 * values are close together or if there is a tiny entry in the Z
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
30 * vector. For each such occurrence the order of the related
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
31 * secular equation problem is reduced by one.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
32 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
33 * DLASD7 is called from DLASD6.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
34 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
35 * Arguments
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
36 * =========
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
37 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
38 * ICOMPQ (input) INTEGER
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
39 * Specifies whether singular vectors are to be computed
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
40 * in compact form, as follows:
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
41 * = 0: Compute singular values only.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
42 * = 1: Compute singular vectors of upper
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
43 * bidiagonal matrix in compact form.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
44 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
45 * NL (input) INTEGER
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
46 * The row dimension of the upper block. NL >= 1.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
47 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
48 * NR (input) INTEGER
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
49 * The row dimension of the lower block. NR >= 1.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
50 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
51 * SQRE (input) INTEGER
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
52 * = 0: the lower block is an NR-by-NR square matrix.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
53 * = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
54 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
55 * The bidiagonal matrix has
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
56 * N = NL + NR + 1 rows and
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
57 * M = N + SQRE >= N columns.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
58 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
59 * K (output) INTEGER
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
60 * Contains the dimension of the non-deflated matrix, this is
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
61 * the order of the related secular equation. 1 <= K <=N.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
62 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
63 * D (input/output) DOUBLE PRECISION array, dimension ( N )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
64 * On entry D contains the singular values of the two submatrices
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
65 * to be combined. On exit D contains the trailing (N-K) updated
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
66 * singular values (those which were deflated) sorted into
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
67 * increasing order.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
68 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
69 * Z (output) DOUBLE PRECISION array, dimension ( M )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
70 * On exit Z contains the updating row vector in the secular
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
71 * equation.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
72 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
73 * ZW (workspace) DOUBLE PRECISION array, dimension ( M )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
74 * Workspace for Z.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
75 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
76 * VF (input/output) DOUBLE PRECISION array, dimension ( M )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
77 * On entry, VF(1:NL+1) contains the first components of all
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
78 * right singular vectors of the upper block; and VF(NL+2:M)
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
79 * contains the first components of all right singular vectors
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
80 * of the lower block. On exit, VF contains the first components
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
81 * of all right singular vectors of the bidiagonal matrix.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
82 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
83 * VFW (workspace) DOUBLE PRECISION array, dimension ( M )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
84 * Workspace for VF.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
85 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
86 * VL (input/output) DOUBLE PRECISION array, dimension ( M )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
87 * On entry, VL(1:NL+1) contains the last components of all
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
88 * right singular vectors of the upper block; and VL(NL+2:M)
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
89 * contains the last components of all right singular vectors
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
90 * of the lower block. On exit, VL contains the last components
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
91 * of all right singular vectors of the bidiagonal matrix.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
92 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
93 * VLW (workspace) DOUBLE PRECISION array, dimension ( M )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
94 * Workspace for VL.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
95 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
96 * ALPHA (input) DOUBLE PRECISION
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
97 * Contains the diagonal element associated with the added row.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
98 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
99 * BETA (input) DOUBLE PRECISION
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
100 * Contains the off-diagonal element associated with the added
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
101 * row.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
102 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
103 * DSIGMA (output) DOUBLE PRECISION array, dimension ( N )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
104 * Contains a copy of the diagonal elements (K-1 singular values
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
105 * and one zero) in the secular equation.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
106 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
107 * IDX (workspace) INTEGER array, dimension ( N )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
108 * This will contain the permutation used to sort the contents of
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
109 * D into ascending order.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
110 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
111 * IDXP (workspace) INTEGER array, dimension ( N )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
112 * This will contain the permutation used to place deflated
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
113 * values of D at the end of the array. On output IDXP(2:K)
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
114 * points to the nondeflated D-values and IDXP(K+1:N)
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
115 * points to the deflated singular values.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
116 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
117 * IDXQ (input) INTEGER array, dimension ( N )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
118 * This contains the permutation which separately sorts the two
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
119 * sub-problems in D into ascending order. Note that entries in
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
120 * the first half of this permutation must first be moved one
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
121 * position backward; and entries in the second half
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
122 * must first have NL+1 added to their values.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
123 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
124 * PERM (output) INTEGER array, dimension ( N )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
125 * The permutations (from deflation and sorting) to be applied
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
126 * to each singular block. Not referenced if ICOMPQ = 0.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
127 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
128 * GIVPTR (output) INTEGER
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
129 * The number of Givens rotations which took place in this
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
130 * subproblem. Not referenced if ICOMPQ = 0.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
131 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
132 * GIVCOL (output) INTEGER array, dimension ( LDGCOL, 2 )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
133 * Each pair of numbers indicates a pair of columns to take place
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
134 * in a Givens rotation. Not referenced if ICOMPQ = 0.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
135 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
136 * LDGCOL (input) INTEGER
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
137 * The leading dimension of GIVCOL, must be at least N.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
138 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
139 * GIVNUM (output) DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
140 * Each number indicates the C or S value to be used in the
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
141 * corresponding Givens rotation. Not referenced if ICOMPQ = 0.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
142 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
143 * LDGNUM (input) INTEGER
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
144 * The leading dimension of GIVNUM, must be at least N.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
145 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
146 * C (output) DOUBLE PRECISION
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
147 * C contains garbage if SQRE =0 and the C-value of a Givens
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
148 * rotation related to the right null space if SQRE = 1.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
149 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
150 * S (output) DOUBLE PRECISION
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
151 * S contains garbage if SQRE =0 and the S-value of a Givens
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
152 * rotation related to the right null space if SQRE = 1.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
153 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
154 * INFO (output) INTEGER
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
155 * = 0: successful exit.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
156 * < 0: if INFO = -i, the i-th argument had an illegal value.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
157 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
158 * Further Details
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
159 * ===============
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
160 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
161 * Based on contributions by
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
162 * Ming Gu and Huan Ren, Computer Science Division, University of
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
163 * California at Berkeley, USA
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
164 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
165 * =====================================================================
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
166 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
167 * .. Parameters ..
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
168 DOUBLE PRECISION ZERO, ONE, TWO, EIGHT
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
169 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0,
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
170 $ EIGHT = 8.0D+0 )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
171 * ..
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
172 * .. Local Scalars ..
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
173 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
174 INTEGER I, IDXI, IDXJ, IDXJP, J, JP, JPREV, K2, M, N,
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
175 $ NLP1, NLP2
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
176 DOUBLE PRECISION EPS, HLFTOL, TAU, TOL, Z1
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
177 * ..
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
178 * .. External Subroutines ..
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
179 EXTERNAL DCOPY, DLAMRG, DROT, XERBLA
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
180 * ..
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
181 * .. External Functions ..
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
182 DOUBLE PRECISION DLAMCH, DLAPY2
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
183 EXTERNAL DLAMCH, DLAPY2
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
184 * ..
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
185 * .. Intrinsic Functions ..
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
186 INTRINSIC ABS, MAX
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
187 * ..
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
188 * .. Executable Statements ..
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
189 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
190 * Test the input parameters.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
191 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
192 INFO = 0
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
193 N = NL + NR + 1
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
194 M = N + SQRE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
195 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
196 IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
197 INFO = -1
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
198 ELSE IF( NL.LT.1 ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
199 INFO = -2
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
200 ELSE IF( NR.LT.1 ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
201 INFO = -3
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
202 ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
203 INFO = -4
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
204 ELSE IF( LDGCOL.LT.N ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
205 INFO = -22
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
206 ELSE IF( LDGNUM.LT.N ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
207 INFO = -24
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
208 END IF
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
209 IF( INFO.NE.0 ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
210 CALL XERBLA( 'DLASD7', -INFO )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
211 RETURN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
212 END IF
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
213 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
214 NLP1 = NL + 1
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
215 NLP2 = NL + 2
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
216 IF( ICOMPQ.EQ.1 ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
217 GIVPTR = 0
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
218 END IF
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
219 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
220 * Generate the first part of the vector Z and move the singular
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
221 * values in the first part of D one position backward.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
222 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
223 Z1 = ALPHA*VL( NLP1 )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
224 VL( NLP1 ) = ZERO
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
225 TAU = VF( NLP1 )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
226 DO 10 I = NL, 1, -1
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
227 Z( I+1 ) = ALPHA*VL( I )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
228 VL( I ) = ZERO
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
229 VF( I+1 ) = VF( I )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
230 D( I+1 ) = D( I )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
231 IDXQ( I+1 ) = IDXQ( I ) + 1
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
232 10 CONTINUE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
233 VF( 1 ) = TAU
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
234 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
235 * Generate the second part of the vector Z.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
236 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
237 DO 20 I = NLP2, M
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
238 Z( I ) = BETA*VF( I )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
239 VF( I ) = ZERO
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
240 20 CONTINUE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
241 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
242 * Sort the singular values into increasing order
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
243 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
244 DO 30 I = NLP2, N
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
245 IDXQ( I ) = IDXQ( I ) + NLP1
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
246 30 CONTINUE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
247 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
248 * DSIGMA, IDXC, IDXC, and ZW are used as storage space.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
249 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
250 DO 40 I = 2, N
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
251 DSIGMA( I ) = D( IDXQ( I ) )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
252 ZW( I ) = Z( IDXQ( I ) )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
253 VFW( I ) = VF( IDXQ( I ) )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
254 VLW( I ) = VL( IDXQ( I ) )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
255 40 CONTINUE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
256 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
257 CALL DLAMRG( NL, NR, DSIGMA( 2 ), 1, 1, IDX( 2 ) )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
258 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
259 DO 50 I = 2, N
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
260 IDXI = 1 + IDX( I )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
261 D( I ) = DSIGMA( IDXI )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
262 Z( I ) = ZW( IDXI )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
263 VF( I ) = VFW( IDXI )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
264 VL( I ) = VLW( IDXI )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
265 50 CONTINUE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
266 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
267 * Calculate the allowable deflation tolerence
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
268 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
269 EPS = DLAMCH( 'Epsilon' )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
270 TOL = MAX( ABS( ALPHA ), ABS( BETA ) )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
271 TOL = EIGHT*EIGHT*EPS*MAX( ABS( D( N ) ), TOL )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
272 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
273 * There are 2 kinds of deflation -- first a value in the z-vector
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
274 * is small, second two (or more) singular values are very close
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
275 * together (their difference is small).
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
276 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
277 * If the value in the z-vector is small, we simply permute the
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
278 * array so that the corresponding singular value is moved to the
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
279 * end.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
280 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
281 * If two values in the D-vector are close, we perform a two-sided
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
282 * rotation designed to make one of the corresponding z-vector
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
283 * entries zero, and then permute the array so that the deflated
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
284 * singular value is moved to the end.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
285 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
286 * If there are multiple singular values then the problem deflates.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
287 * Here the number of equal singular values are found. As each equal
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
288 * singular value is found, an elementary reflector is computed to
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
289 * rotate the corresponding singular subspace so that the
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
290 * corresponding components of Z are zero in this new basis.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
291 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
292 K = 1
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
293 K2 = N + 1
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
294 DO 60 J = 2, N
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
295 IF( ABS( Z( J ) ).LE.TOL ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
296 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
297 * Deflate due to small z component.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
298 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
299 K2 = K2 - 1
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
300 IDXP( K2 ) = J
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
301 IF( J.EQ.N )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
302 $ GO TO 100
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
303 ELSE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
304 JPREV = J
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
305 GO TO 70
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
306 END IF
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
307 60 CONTINUE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
308 70 CONTINUE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
309 J = JPREV
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
310 80 CONTINUE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
311 J = J + 1
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
312 IF( J.GT.N )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
313 $ GO TO 90
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
314 IF( ABS( Z( J ) ).LE.TOL ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
315 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
316 * Deflate due to small z component.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
317 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
318 K2 = K2 - 1
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
319 IDXP( K2 ) = J
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
320 ELSE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
321 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
322 * Check if singular values are close enough to allow deflation.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
323 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
324 IF( ABS( D( J )-D( JPREV ) ).LE.TOL ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
325 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
326 * Deflation is possible.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
327 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
328 S = Z( JPREV )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
329 C = Z( J )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
330 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
331 * Find sqrt(a**2+b**2) without overflow or
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
332 * destructive underflow.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
333 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
334 TAU = DLAPY2( C, S )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
335 Z( J ) = TAU
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
336 Z( JPREV ) = ZERO
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
337 C = C / TAU
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
338 S = -S / TAU
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
339 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
340 * Record the appropriate Givens rotation
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
341 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
342 IF( ICOMPQ.EQ.1 ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
343 GIVPTR = GIVPTR + 1
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
344 IDXJP = IDXQ( IDX( JPREV )+1 )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
345 IDXJ = IDXQ( IDX( J )+1 )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
346 IF( IDXJP.LE.NLP1 ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
347 IDXJP = IDXJP - 1
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
348 END IF
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
349 IF( IDXJ.LE.NLP1 ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
350 IDXJ = IDXJ - 1
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
351 END IF
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
352 GIVCOL( GIVPTR, 2 ) = IDXJP
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
353 GIVCOL( GIVPTR, 1 ) = IDXJ
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
354 GIVNUM( GIVPTR, 2 ) = C
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
355 GIVNUM( GIVPTR, 1 ) = S
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
356 END IF
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
357 CALL DROT( 1, VF( JPREV ), 1, VF( J ), 1, C, S )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
358 CALL DROT( 1, VL( JPREV ), 1, VL( J ), 1, C, S )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
359 K2 = K2 - 1
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
360 IDXP( K2 ) = JPREV
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
361 JPREV = J
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
362 ELSE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
363 K = K + 1
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
364 ZW( K ) = Z( JPREV )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
365 DSIGMA( K ) = D( JPREV )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
366 IDXP( K ) = JPREV
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
367 JPREV = J
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
368 END IF
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
369 END IF
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
370 GO TO 80
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
371 90 CONTINUE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
372 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
373 * Record the last singular value.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
374 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
375 K = K + 1
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
376 ZW( K ) = Z( JPREV )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
377 DSIGMA( K ) = D( JPREV )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
378 IDXP( K ) = JPREV
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
379 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
380 100 CONTINUE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
381 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
382 * Sort the singular values into DSIGMA. The singular values which
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
383 * were not deflated go into the first K slots of DSIGMA, except
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
384 * that DSIGMA(1) is treated separately.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
385 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
386 DO 110 J = 2, N
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
387 JP = IDXP( J )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
388 DSIGMA( J ) = D( JP )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
389 VFW( J ) = VF( JP )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
390 VLW( J ) = VL( JP )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
391 110 CONTINUE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
392 IF( ICOMPQ.EQ.1 ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
393 DO 120 J = 2, N
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
394 JP = IDXP( J )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
395 PERM( J ) = IDXQ( IDX( JP )+1 )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
396 IF( PERM( J ).LE.NLP1 ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
397 PERM( J ) = PERM( J ) - 1
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
398 END IF
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
399 120 CONTINUE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
400 END IF
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
401 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
402 * The deflated singular values go back into the last N - K slots of
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
403 * D.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
404 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
405 CALL DCOPY( N-K, DSIGMA( K+1 ), 1, D( K+1 ), 1 )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
406 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
407 * Determine DSIGMA(1), DSIGMA(2), Z(1), VF(1), VL(1), VF(M), and
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
408 * VL(M).
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
409 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
410 DSIGMA( 1 ) = ZERO
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
411 HLFTOL = TOL / TWO
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
412 IF( ABS( DSIGMA( 2 ) ).LE.HLFTOL )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
413 $ DSIGMA( 2 ) = HLFTOL
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
414 IF( M.GT.N ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
415 Z( 1 ) = DLAPY2( Z1, Z( M ) )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
416 IF( Z( 1 ).LE.TOL ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
417 C = ONE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
418 S = ZERO
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
419 Z( 1 ) = TOL
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
420 ELSE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
421 C = Z1 / Z( 1 )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
422 S = -Z( M ) / Z( 1 )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
423 END IF
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
424 CALL DROT( 1, VF( M ), 1, VF( 1 ), 1, C, S )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
425 CALL DROT( 1, VL( M ), 1, VL( 1 ), 1, C, S )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
426 ELSE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
427 IF( ABS( Z1 ).LE.TOL ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
428 Z( 1 ) = TOL
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
429 ELSE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
430 Z( 1 ) = Z1
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
431 END IF
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
432 END IF
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
433 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
434 * Restore Z, VF, and VL.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
435 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
436 CALL DCOPY( K-1, ZW( 2 ), 1, Z( 2 ), 1 )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
437 CALL DCOPY( N-1, VFW( 2 ), 1, VF( 2 ), 1 )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
438 CALL DCOPY( N-1, VLW( 2 ), 1, VL( 2 ), 1 )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
439 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
440 RETURN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
441 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
442 * End of DLASD7
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
443 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
444 END