Mercurial > octave-nkf
annotate libcruft/lapack/slarft.f @ 7948:af10baa63915 ss-3-1-50
3.1.50 snapshot
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Fri, 18 Jul 2008 17:42:48 -0400 |
parents | 82be108cc558 |
children |
rev | line source |
---|---|
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
1 SUBROUTINE SLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT ) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
2 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
3 * -- LAPACK auxiliary routine (version 3.1) -- |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
4 * 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
|
5 * November 2006 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
6 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
7 * .. Scalar Arguments .. |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
8 CHARACTER DIRECT, STOREV |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
9 INTEGER K, LDT, LDV, N |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
10 * .. |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
11 * .. Array Arguments .. |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
12 REAL T( LDT, * ), TAU( * ), V( LDV, * ) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
13 * .. |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
14 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
15 * Purpose |
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 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
18 * SLARFT forms the triangular factor T of a real block reflector H |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
19 * of order n, which is defined as a product of k elementary reflectors. |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
20 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
21 * If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
22 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
23 * If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular. |
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 * If STOREV = 'C', the vector which defines the elementary reflector |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
26 * H(i) is stored in the i-th column of the array V, and |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
27 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
28 * H = I - V * T * V' |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
29 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
30 * If STOREV = 'R', the vector which defines the elementary reflector |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
31 * H(i) is stored in the i-th row of the array V, and |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
32 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
33 * H = I - V' * T * V |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
34 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
35 * Arguments |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
36 * ========= |
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 * DIRECT (input) CHARACTER*1 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
39 * Specifies the order in which the elementary reflectors are |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
40 * multiplied to form the block reflector: |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
41 * = 'F': H = H(1) H(2) . . . H(k) (Forward) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
42 * = 'B': H = H(k) . . . H(2) H(1) (Backward) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
43 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
44 * STOREV (input) CHARACTER*1 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
45 * Specifies how the vectors which define the elementary |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
46 * reflectors are stored (see also Further Details): |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
47 * = 'C': columnwise |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
48 * = 'R': rowwise |
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 * N (input) INTEGER |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
51 * The order of the block reflector H. N >= 0. |
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 * K (input) INTEGER |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
54 * The order of the triangular factor T (= the number of |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
55 * elementary reflectors). K >= 1. |
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 * V (input/output) REAL array, dimension |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
58 * (LDV,K) if STOREV = 'C' |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
59 * (LDV,N) if STOREV = 'R' |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
60 * The matrix V. See further details. |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
61 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
62 * LDV (input) INTEGER |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
63 * The leading dimension of the array V. |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
64 * If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K. |
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 * TAU (input) REAL array, dimension (K) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
67 * TAU(i) must contain the scalar factor of the elementary |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
68 * reflector H(i). |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
69 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
70 * T (output) REAL array, dimension (LDT,K) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
71 * The k by k triangular factor T of the block reflector. |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
72 * If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
73 * lower triangular. The rest of the array is not used. |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
74 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
75 * LDT (input) INTEGER |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
76 * The leading dimension of the array T. LDT >= K. |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
77 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
78 * Further Details |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
79 * =============== |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
80 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
81 * The shape of the matrix V and the storage of the vectors which define |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
82 * the H(i) is best illustrated by the following example with n = 5 and |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
83 * k = 3. The elements equal to 1 are not stored; the corresponding |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
84 * array elements are modified but restored on exit. The rest of the |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
85 * array is not used. |
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 * DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and STOREV = 'R': |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
88 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
89 * V = ( 1 ) V = ( 1 v1 v1 v1 v1 ) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
90 * ( v1 1 ) ( 1 v2 v2 v2 ) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
91 * ( v1 v2 1 ) ( 1 v3 v3 ) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
92 * ( v1 v2 v3 ) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
93 * ( v1 v2 v3 ) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
94 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
95 * DIRECT = 'B' and STOREV = 'C': DIRECT = 'B' and STOREV = 'R': |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
96 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
97 * V = ( v1 v2 v3 ) V = ( v1 v1 1 ) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
98 * ( v1 v2 v3 ) ( v2 v2 v2 1 ) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
99 * ( 1 v2 v3 ) ( v3 v3 v3 v3 1 ) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
100 * ( 1 v3 ) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
101 * ( 1 ) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
102 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
103 * ===================================================================== |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
104 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
105 * .. Parameters .. |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
106 REAL ONE, ZERO |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
107 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
108 * .. |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
109 * .. Local Scalars .. |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
110 INTEGER I, J |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
111 REAL VII |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
112 * .. |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
113 * .. External Subroutines .. |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
114 EXTERNAL SGEMV, STRMV |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
115 * .. |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
116 * .. External Functions .. |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
117 LOGICAL LSAME |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
118 EXTERNAL LSAME |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
119 * .. |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
120 * .. Executable Statements .. |
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 * Quick return if possible |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
123 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
124 IF( N.EQ.0 ) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
125 $ RETURN |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
126 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
127 IF( LSAME( DIRECT, 'F' ) ) THEN |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
128 DO 20 I = 1, K |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
129 IF( TAU( I ).EQ.ZERO ) THEN |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
130 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
131 * H(i) = I |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
132 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
133 DO 10 J = 1, I |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
134 T( J, I ) = ZERO |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
135 10 CONTINUE |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
136 ELSE |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
137 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
138 * general case |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
139 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
140 VII = V( I, I ) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
141 V( I, I ) = ONE |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
142 IF( LSAME( STOREV, 'C' ) ) THEN |
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 * T(1:i-1,i) := - tau(i) * V(i:n,1:i-1)' * V(i:n,i) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
145 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
146 CALL SGEMV( 'Transpose', N-I+1, I-1, -TAU( I ), |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
147 $ V( I, 1 ), LDV, V( I, I ), 1, ZERO, |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
148 $ T( 1, I ), 1 ) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
149 ELSE |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
150 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
151 * T(1:i-1,i) := - tau(i) * V(1:i-1,i:n) * V(i,i:n)' |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
152 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
153 CALL SGEMV( 'No transpose', I-1, N-I+1, -TAU( I ), |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
154 $ V( 1, I ), LDV, V( I, I ), LDV, ZERO, |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
155 $ T( 1, I ), 1 ) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
156 END IF |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
157 V( I, I ) = VII |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
158 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
159 * T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
160 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
161 CALL STRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T, |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
162 $ LDT, T( 1, I ), 1 ) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
163 T( I, I ) = TAU( I ) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
164 END IF |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
165 20 CONTINUE |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
166 ELSE |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
167 DO 40 I = K, 1, -1 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
168 IF( TAU( I ).EQ.ZERO ) THEN |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
169 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
170 * H(i) = I |
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 DO 30 J = I, K |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
173 T( J, I ) = ZERO |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
174 30 CONTINUE |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
175 ELSE |
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 * general case |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
178 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
179 IF( I.LT.K ) THEN |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
180 IF( LSAME( STOREV, 'C' ) ) THEN |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
181 VII = V( N-K+I, I ) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
182 V( N-K+I, I ) = ONE |
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 * T(i+1:k,i) := |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
185 * - tau(i) * V(1:n-k+i,i+1:k)' * V(1:n-k+i,i) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
186 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
187 CALL SGEMV( 'Transpose', N-K+I, K-I, -TAU( I ), |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
188 $ V( 1, I+1 ), LDV, V( 1, I ), 1, ZERO, |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
189 $ T( I+1, I ), 1 ) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
190 V( N-K+I, I ) = VII |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
191 ELSE |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
192 VII = V( I, N-K+I ) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
193 V( I, N-K+I ) = ONE |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
194 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
195 * T(i+1:k,i) := |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
196 * - tau(i) * V(i+1:k,1:n-k+i) * V(i,1:n-k+i)' |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
197 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
198 CALL SGEMV( 'No transpose', K-I, N-K+I, -TAU( I ), |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
199 $ V( I+1, 1 ), LDV, V( I, 1 ), LDV, ZERO, |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
200 $ T( I+1, I ), 1 ) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
201 V( I, N-K+I ) = VII |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
202 END IF |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
203 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
204 * T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
205 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
206 CALL STRMV( 'Lower', 'No transpose', 'Non-unit', K-I, |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
207 $ T( I+1, I+1 ), LDT, T( I+1, I ), 1 ) |
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 T( I, I ) = TAU( I ) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
210 END IF |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
211 40 CONTINUE |
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 RETURN |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
214 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
215 * End of SLARFT |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
216 * |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff
changeset
|
217 END |