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
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
|
SUBROUTINE EREDUC(E, M, N, Q, Z, ISTAIR, RANKE, TOL)
C PURPOSE:
C
C Given an M x N matrix E (not necessarily regular) the subroutine
C EREDUC computes a unitary transformed matrix Q*E*Z which is in
C column echelon form (trapezoidal form). Furthermore the rank of
C matrix E is determined.
C
C CONTRIBUTOR: Th.G.J. Beelen (Philips Glass Eindhoven).
C Copyright SLICOT
C
C REVISIONS: 1988, January 29.
C
C Specification of the parameters.
C
C .. Scalar arguments ..
C
INTEGER LDE, LDQ, LDZ, M, N, RANKE
DOUBLE PRECISION TOL
C
C .. Array arguments ..
C
INTEGER ISTAIR(M)
C DOUBLE PRECISION E(LDE,N), Q(LDQ,M), Z(LDZ,N)
C SET E(M,N) Q(M,M) Z(N,N)
DOUBLE PRECISION E(M,N), Q(M,M), Z(N,N)
C Local variables.
C
INTEGER I, J, JMX, K, KM1, L, LK, MNK, NR1
DOUBLE PRECISION EMXNRM, EMX, SC, SS
LOGICAL LZERO
C
LDE=M
LDQ=M
LDZ=N
do 991 i=1,m
do 991 j=1,m
q(i,j)=0.d0
991 continue
do 992 i=1,m
q(i,i)=1.0d0
992 continue
do 993 i=1,n
do 993 j=1,n
z(i,j)=0.d0
993 continue
do 994 i=1,n
z(i,i)=1.0d0
994 continue
RANKE = MIN0(M,N)
C
K = N
LZERO = .FALSE.
C
C WHILE ((K > 0) AND (NOT a zero submatrix encountered) DO
10 IF ((K .GT. 0) .AND. (.NOT.LZERO)) THEN
C
C
MNK = M - N + K
EMXNRM = 0.0D0
LK = MNK
DO 20 L = MNK, 1, -1
JMX = IDAMAX(K, E(L,1), LDE)
EMX = DABS(E(L,JMX))
IF (EMX .GT. EMXNRM) THEN
EMXNRM = EMX
LK = L
END IF
20 CONTINUE
C
IF (EMXNRM .LT. TOL) THEN
C
C Set submatrix Ek to zero.
C
DO 40 J = 1, K
DO 30 L = 1, MNK
E(L,J) = 0.0D0
30 CONTINUE
40 CONTINUE
LZERO = .TRUE.
RANKE = N - K
ELSE
C
C Submatrix Ek is not considered to be identically zero.
C Check whether rows have to be interchanged.
C
IF (LK .NE. MNK) THEN
C
C Interchange rows lk and m-n+k in whole E-matrix and
C update the row transformation matrix Q.
C (# elements involved = m)
C
CALL DSWAP(N, E(LK,1), LDE, E(MNK,1), LDE)
CALL DSWAP(M, Q(LK,1), LDQ, Q(MNK,1), LDQ)
END IF
C
KM1 = K - 1
DO 50 J = 1, KM1
C
C Determine the column Givens transformation to annihilate
C E(m-n+k,j) using E(m-n+k,k) as pivot.
C Apply the transformation to the columns of Ek.
C (# elements involved = m-n+k)
C Update the column transformation matrix Z.
C (# elements involved = n)
C
CALL DGIV(E(MNK,K), E(MNK,J), SC, SS)
CALL DROT(MNK, E(1,K), 1, E(1,J), 1, SC, SS)
E(MNK, J) = 0.0D0
CALL DROT(N, Z(1,K), 1, Z(1,J), 1, SC, SS)
50 CONTINUE
C
K = K - 1
END IF
GOTO 10
END IF
C END WHILE 10
C
C Initialise administration staircase form, i.e.,
C ISTAIR(i) = j if E(i,j) is a nonzero corner point
C = -j if E(i,j) is on the boundary but is no corner pt.
C Thus,
C ISTAIR(m-k) = n-k for k=0,...,rank(E)-1
C = -(n-rank(E)+1) for k=rank(E),...,m-1.
C
DO 60 I = 1, RANKE
ISTAIR(M - I + 1) = N - I + 1
60 CONTINUE
C
NR1 = N - RANKE + 1
DO 70 I = RANKE, M - 1
ISTAIR(M-I) = -NR1
70 CONTINUE
C
RETURN
C *** Last line of EREDUC *********************************************
END
|
