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
|
SUBROUTINE R1MPYQ(M,N,A,LDA,V,W)
INTEGER M,N,LDA
DOUBLE PRECISION A(LDA,N),V(N),W(N)
C **********
C
C SUBROUTINE R1MPYQ
C
C GIVEN AN M BY N MATRIX A, THIS SUBROUTINE COMPUTES A*Q WHERE
C Q IS THE PRODUCT OF 2*(N - 1) TRANSFORMATIONS
C
C GV(N-1)*...*GV(1)*GW(1)*...*GW(N-1)
C
C AND GV(I), GW(I) ARE GIVENS ROTATIONS IN THE (I,N) PLANE WHICH
C ELIMINATE ELEMENTS IN THE I-TH AND N-TH PLANES, RESPECTIVELY.
C Q ITSELF IS NOT GIVEN, RATHER THE INFORMATION TO RECOVER THE
C GV, GW ROTATIONS IS SUPPLIED.
C
C THE SUBROUTINE STATEMENT IS
C
C SUBROUTINE R1MPYQ(M,N,A,LDA,V,W)
C
C WHERE
C
C M IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
C OF ROWS OF A.
C
C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
C OF COLUMNS OF A.
C
C A IS AN M BY N ARRAY. ON INPUT A MUST CONTAIN THE MATRIX
C TO BE POSTMULTIPLIED BY THE ORTHOGONAL MATRIX Q
C DESCRIBED ABOVE. ON OUTPUT A*Q HAS REPLACED A.
C
C LDA IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN M
C WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY A.
C
C V IS AN INPUT ARRAY OF LENGTH N. V(I) MUST CONTAIN THE
C INFORMATION NECESSARY TO RECOVER THE GIVENS ROTATION GV(I)
C DESCRIBED ABOVE.
C
C W IS AN INPUT ARRAY OF LENGTH N. W(I) MUST CONTAIN THE
C INFORMATION NECESSARY TO RECOVER THE GIVENS ROTATION GW(I)
C DESCRIBED ABOVE.
C
C SUBROUTINES CALLED
C
C FORTRAN-SUPPLIED ... DABS,DSQRT
C
C MINPACK. VERSION OF DECEMBER 1978.
C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE
C
C **********
INTEGER I,J,NMJ,NM1
DOUBLE PRECISION COS,ONE,SIN,TEMP
DATA ONE /1.0D0/
C
C APPLY THE FIRST SET OF GIVENS ROTATIONS TO A.
C
NM1 = N - 1
IF (NM1 .LT. 1) GO TO 50
DO 20 NMJ = 1, NM1
J = N - NMJ
IF (DABS(V(J)) .GT. ONE) COS = ONE/V(J)
IF (DABS(V(J)) .GT. ONE) SIN = DSQRT(ONE-COS**2)
IF (DABS(V(J)) .LE. ONE) SIN = V(J)
IF (DABS(V(J)) .LE. ONE) COS = DSQRT(ONE-SIN**2)
DO 10 I = 1, M
TEMP = COS*A(I,J) - SIN*A(I,N)
A(I,N) = SIN*A(I,J) + COS*A(I,N)
A(I,J) = TEMP
10 CONTINUE
20 CONTINUE
C
C APPLY THE SECOND SET OF GIVENS ROTATIONS TO A.
C
DO 40 J = 1, NM1
IF (DABS(W(J)) .GT. ONE) COS = ONE/W(J)
IF (DABS(W(J)) .GT. ONE) SIN = DSQRT(ONE-COS**2)
IF (DABS(W(J)) .LE. ONE) SIN = W(J)
IF (DABS(W(J)) .LE. ONE) COS = DSQRT(ONE-SIN**2)
DO 30 I = 1, M
TEMP = COS*A(I,J) + SIN*A(I,N)
A(I,N) = -SIN*A(I,J) + COS*A(I,N)
A(I,J) = TEMP
30 CONTINUE
40 CONTINUE
50 CONTINUE
RETURN
C
C LAST CARD OF SUBROUTINE R1MPYQ.
C
END
|
