1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
|
*DECK ALH12
SUBROUTINE ALH12(MODE,LPIVOT,L1,M,U,UP,C,ICE,ICV,NCV)
*
*-----------------------------------------------------------------------
*
*Purpose:
* construction and/or application of a single Householder
* transformation (Q = I + U*(U**T)/B).
*
*Author(s): C. L. Lawson and R. J. Hanson
*
*Parameters: input
* MODE algorithm flag (=1/2: Selects algorithm H1 to construct and
* apply a Householder transformation / Algorithm H2 to apply a
* previously constructed transformation).
* LPIVOT index of the pivot element.
* L1,M if L1.le.M, the transformation will be constructed to
* zero elements indexed from L1 through M. If L1.gt.M,
* the subroutine does an identity transformation.
* U pivot vector (if MODE=1). At exit, contains quantities
* defining the vector U of the Householder transformation.
* On entry with MODE=2, U should contain information
* previously computed with MODE=1.
* C matrix which will be regarded as a set of vectors to which the
* Householder transformation is to be applied.
* ICE storage increment between elements of vectors in C.
* ICV storage increment between vectors in C.
* NCV number of vectors in C to be transformed. if NCV.le.0, no
* operations will be done on C.
*
*Parameters: output
* U quantitie defining the vector U of the Householder
* transformation. These will not be modified during the entry
* with MODE=2.
* UP quantity defining the vector U of the Householder
* transformation. On entry with MODE = 2, UP should
* contain information previously computed with MODE=1. These
* will not be modified during the entry with MODE=2.
* C set of transformed vectors.
*
*-----------------------------------------------------------------------
*
IMPLICIT NONE
*----
* Subroutine arguments
*----
INTEGER MODE, LPIVOT, L1, M, ICE, ICV, NCV
DOUBLE PRECISION U(M),C(M),UP
*----
* Local variables
*----
INTEGER I, I2, I3, I4, INCR, J
DOUBLE PRECISION B, CL, CLINV, SM
*
IF(0.GE.LPIVOT .OR. LPIVOT.GE.L1 .OR. L1.GT.M) RETURN
CL = ABS(U(LPIVOT))
IF(MODE.NE.2) THEN
*----
* Construct the transformation.
*----
DO J = L1, M
CL=MAX(ABS(U(J)),CL)
ENDDO
IF(CL.LE.0) RETURN
CLINV = 1.0D0 / CL
SM = (U(LPIVOT)*CLINV) ** 2 + SUM( (U(L1:M)*CLINV)**2 )
CL = CL * SQRT(SM)
IF(U(LPIVOT).GT.0) THEN
CL = -CL
ENDIF
UP = U(LPIVOT) - CL
U(LPIVOT) = CL
ELSE
*----
* Apply the transformation I+U*(U**T)/B to C.
*----
IF(CL.LE.0) RETURN
ENDIF
IF(NCV.LE.0) RETURN
B = UP * U(LPIVOT)
*----
* B must be nonpositive here. If B = 0., return.
*----
IF(B.LT.0.0) THEN
B = 1.0D0 / B
I2 = 1 - ICV + ICE * (LPIVOT-1)
INCR = ICE * (L1-LPIVOT)
DO J = 1, NCV
I2 = I2 + ICV
I3 = I2 + INCR
I4 = I3
SM = C(I2) * UP
DO I = L1, M
SM = SM + C(I3) * U(I)
I3 = I3 + ICE
ENDDO
IF(SM.NE.0) THEN
SM = SM * B
C(I2) = C(I2) + SM * UP
DO I = L1, M
C(I4) = C(I4) + SM * U(I)
I4 = I4 + ICE
ENDDO
ENDIF
ENDDO
ENDIF
RETURN
END
|
