]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 1 | * |
2 | * $Id$ | |
3 | * | |
4 | * $Log$ | |
5 | * Revision 1.1.1.1 1996/04/01 15:02:27 mclareni | |
6 | * Mathlib gen | |
7 | * | |
8 | * | |
9 | #include "gen/pilot.h" | |
10 | SUBROUTINE CHEB(M,N,MDIM,NDIM,A,B,TOL,RELERR,X,RANK,RESMAX, | |
11 | + ITER,OCODE) | |
12 | C LINEAR CHEBYSHEV FIT ROUTINE FROM ACM TOMS 1, 266(1975) | |
13 | C BY BARRODALE AND PHILLIPS | |
14 | DIMENSION A(NDIM,MDIM), B(MDIM), X(NDIM) | |
15 | INTEGER PROW,PCOL,RANK,RANKP1,OCODE | |
16 | #include "e221prec.inc" | |
17 | MP1 = M + 1 | |
18 | NP1 = N + 1 | |
19 | NP2 = N + 2 | |
20 | NP3 = N + 3 | |
21 | NP1MR = 1 | |
22 | RANK = N | |
23 | RELTMP = RELERR | |
24 | RELERR = 0. | |
25 | DO 10 J= 1,M | |
26 | A(NP1,J) = 1. | |
27 | A(NP2,J) = -B(J) | |
28 | A(NP3,J) = N+J | |
29 | 10 CONTINUE | |
30 | A(NP1,MP1) = 0. | |
31 | ITER = 0 | |
32 | OCODE = 1 | |
33 | DO 20 I= 1, N | |
34 | X(I)=0. | |
35 | A(I,MP1)=I | |
36 | 20 CONTINUE | |
37 | C LEVEL 1. | |
38 | LEV = 1 | |
39 | K = 0 | |
40 | 30 K = K + 1 | |
41 | KP1 = K + 1 | |
42 | NP1MK = NP1 - K | |
43 | MODE = 0 | |
44 | CALL VFILL(B(K),M-K+1,1.) | |
45 | C DETERMINE THE VECTOR TO ENTER THE BASIS. | |
46 | 50 D = -BIG | |
47 | DO 60 J= K, M | |
48 | IF (B(J) .EQ. 0.) GOTO 60 | |
49 | DD = ABS(A(NP2,J)) | |
50 | IF (DD .LE. D) GOTO 60 | |
51 | PCOL = J | |
52 | D = DD | |
53 | 60 CONTINUE | |
54 | IF (K .GT. 1) GOTO 70 | |
55 | C TEST FOR ZERO RIGHT-HAND SIDE | |
56 | IF (D .GT. TOL) GOTO 70 | |
57 | RESMAX = 0. | |
58 | MODE = 2 | |
59 | GOTO 380 | |
60 | C DETERMINE THE VECTOR TO LEAVE THE BASIS. | |
61 | 70 D = TOL | |
62 | DO 80 I= 1, NP1MK | |
63 | DD = ABS(A(I,PCOL)) | |
64 | IF (DD .LE. D) GOTO 80 | |
65 | PROW = I | |
66 | D = DD | |
67 | 80 CONTINUE | |
68 | IF (D .GT. TOL) GOTO 330 | |
69 | C CHECK FOR LINEAR DEPENDENCE IN LEVEL 1. | |
70 | B(PCOL) = 0. | |
71 | IF (MODE .EQ. 1) GOTO 50 | |
72 | DO 100 J= K, M | |
73 | IF (B(J) .EQ. 0.) GOTO 100 | |
74 | DO 90 I= 1, NP1MK | |
75 | IF (ABS(A(I,J)) .LE. TOL) GOTO 90 | |
76 | MODE = 1 | |
77 | GOTO 50 | |
78 | 90 CONTINUE | |
79 | 100 CONTINUE | |
80 | RANK = K - 1 | |
81 | NP1MR = NP1 - RANK | |
82 | OCODE = 0 | |
83 | GOTO 160 | |
84 | 110 IF (PCOL .EQ. K) GOTO 130 | |
85 | C INTERCHANGE COLUMNS IN LEVEL 1. | |
86 | DO 120 I= 1, NP3 | |
87 | D = A(I,PCOL) | |
88 | A(I,PCOL) = A(I,K) | |
89 | A(I,K) = D | |
90 | 120 CONTINUE | |
91 | 130 IF (PROW .EQ. NP1MK) GOTO 150 | |
92 | C INTERCHANGE ROWS IN LEVEL 1 | |
93 | DO 140 J= 1, MP1 | |
94 | D = A(PROW,J) | |
95 | A(PROW,J) = A(NP1MK,J) | |
96 | A(NP1MK,J) = D | |
97 | 140 CONTINUE | |
98 | 150 IF (K.LT.N) GOTO 30 | |
99 | 160 IF (RANK.EQ.M) GOTO 380 | |
100 | RANKP1 = RANK + 1 | |
101 | C LEVEL 2. | |
102 | LEV = 2 | |
103 | C DETERMINE THE VECTOR TO ENTER THE BASIS. | |
104 | D = TOL | |
105 | DO 170 J= RANKP1, M | |
106 | DD = ABS(A(NP2,J)) | |
107 | IF (DD .LE. D) GOTO 170 | |
108 | PCOL = J | |
109 | D = DD | |
110 | 170 CONTINUE | |
111 | C COMPARE CHEBYSHEV ERROR WITH TOL. | |
112 | IF (D.GT.TOL) GOTO 180 | |
113 | RESMAX = 0. | |
114 | MODE = 3 | |
115 | GOTO 380 | |
116 | C | |
117 | 180 IF (A(NP2,PCOL) .LT. -TOL) GOTO 200 | |
118 | A(NP1,PCOL) = 2. - A(NP1,PCOL) | |
119 | DO 190 I= NP1MR, NP3 | |
120 | IF (I.EQ.NP1) GOTO 190 | |
121 | A(I,PCOL) = -A(I,PCOL) | |
122 | 190 CONTINUE | |
123 | C ARRANGE FOR ALL ENTRIES IN PIVOT COL (EXC. PIVOT) TO BE NEGATIVE | |
124 | 200 DO 220 I= NP1MR, N | |
125 | IF (A(I,PCOL) .LT. TOL) GOTO 220 | |
126 | DO 210 J= 1, M | |
127 | A(NP1,J) = A(NP1,J) + 2.*A(I,J) | |
128 | A(I,J) = -A(I,J) | |
129 | 210 CONTINUE | |
130 | A(I,MP1) = -A(I,MP1) | |
131 | 220 CONTINUE | |
132 | PROW = NP1 | |
133 | GOTO 330 | |
134 | C | |
135 | 230 IF (RANKP1 .EQ. M) GOTO 380 | |
136 | IF (PCOL .EQ. M) GOTO 250 | |
137 | C INTERCHANGE COLUMNS IN LEVEL 2. | |
138 | DO 240 I= NP1MR,NP3 | |
139 | D = A(I,PCOL) | |
140 | A(I,PCOL) =A(I,M) | |
141 | A(I,M) = D | |
142 | 240 CONTINUE | |
143 | 250 MM1 = M-1 | |
144 | C LEVEL 3. | |
145 | LEV = 3 | |
146 | C DETERMINE THE VECTOR TO ENTER THE BASIS. | |
147 | 260 D= -TOL | |
148 | VAL = 2. * A(NP2,M) | |
149 | DO 280 J= RANKP1,MM1 | |
150 | IF (A(NP2,J) .GE. D) GOTO 270 | |
151 | PCOL = J | |
152 | D = A(NP2,J) | |
153 | MODE = 0 | |
154 | GOTO 280 | |
155 | 270 DD = VAL - A(NP2,J) | |
156 | IF (DD .GE. D) GOTO 280 | |
157 | MODE = 1 | |
158 | PCOL = J | |
159 | D = DD | |
160 | 280 CONTINUE | |
161 | IF (D .GE. -TOL) GOTO 380 | |
162 | DD = -D/A(NP2,M) | |
163 | IF (DD .GE. RELTMP) GOTO 290 | |
164 | RELERR = DD | |
165 | MODE = 4 | |
166 | GOTO 380 | |
167 | C | |
168 | 290 IF (MODE .EQ. 0) GOTO 310 | |
169 | DO 300 I= NP1MR,NP1 | |
170 | A(I,PCOL) = 2.*A(I,M) - A(I,PCOL) | |
171 | 300 CONTINUE | |
172 | A(NP2,PCOL) = D | |
173 | A(NP3,PCOL) = -A(NP3,PCOL) | |
174 | C DETERMINE THE VECTOR TO LEAVE THE BASIS. | |
175 | 310 D = BIG | |
176 | DO 320 I= NP1MR,NP1 | |
177 | IF (A(I,PCOL) .LE. TOL) GOTO 320 | |
178 | DD = A(I,M)/A(I,PCOL) | |
179 | IF (DD .GE. D) GOTO 320 | |
180 | PROW = I | |
181 | D = DD | |
182 | 320 CONTINUE | |
183 | IF (D .LT. BIG) GOTO 330 | |
184 | OCODE = 2 | |
185 | GOTO 380 | |
186 | C PIVOT ON A(PROW,PCOL) | |
187 | 330 PIVOT = A(PROW,PCOL) | |
188 | DO 340 J= 1, M | |
189 | A(PROW,J) = A(PROW,J)/PIVOT | |
190 | 340 CONTINUE | |
191 | DO 360 J= 1, M | |
192 | IF (J.EQ.PCOL) GOTO 360 | |
193 | D = A(PROW,J) | |
194 | DO 350 I= NP1MR,NP2 | |
195 | IF (I.EQ.PROW) GOTO 350 | |
196 | A(I,J) = A(I,J) - D*A(I,PCOL) | |
197 | 350 CONTINUE | |
198 | 360 CONTINUE | |
199 | TPIVOT = -PIVOT | |
200 | DO 370 I= NP1MR, NP2 | |
201 | A(I,PCOL) = A(I,PCOL)/TPIVOT | |
202 | 370 CONTINUE | |
203 | A(PROW,PCOL) = 1./PIVOT | |
204 | D = A(PROW,MP1) | |
205 | A(PROW,MP1) = A(NP3,PCOL) | |
206 | A(NP3,PCOL) = D | |
207 | ITER = ITER + 1 | |
208 | GOTO (110, 230, 260), LEV | |
209 | C | |
210 | C PREPARE OUTPUT | |
211 | 380 DO 390 J= 1, M | |
212 | B(J) = 0. | |
213 | 390 CONTINUE | |
214 | IF (MODE .EQ. 2) GOTO 450 | |
215 | DO 400 J= 1, RANK | |
216 | K= A(NP3,J) | |
217 | X(K) = A(NP2,J) | |
218 | 400 CONTINUE | |
219 | IF (MODE.EQ.3 .OR. RANK.EQ.M) GOTO 450 | |
220 | DO 410 I= NP1MR,NP1 | |
221 | K = ABS(A(I,MP1)) - N | |
222 | B(K) = A(NP2,M) * SIGN(1., A(I,MP1)) | |
223 | 410 CONTINUE | |
224 | IF (RANKP1 .EQ. M) GOTO 430 | |
225 | DO 420 J= RANKP1, MM1 | |
226 | K = ABS(A(NP3,J)) - N | |
227 | B(K) = (A(NP2,M)-A(NP2,J)) * SIGN(1., A(NP3,J)) | |
228 | 420 CONTINUE | |
229 | C TEST FOR NON-UNIQUE SOLUTION | |
230 | 430 DO 440 I= NP1MR, NP1 | |
231 | IF (ABS(A(I,M)) .GT. TOL) GOTO 440 | |
232 | OCODE = 0 | |
233 | GOTO 450 | |
234 | 440 CONTINUE | |
235 | 450 IF (MODE.NE.2 .AND. MODE.NE.3) RESMAX = A(NP2,M) | |
236 | IF (RANK.EQ.M) RESMAX=0. | |
237 | IF (MODE.EQ.4) RESMAX = RESMAX-D | |
238 | RETURN | |
239 | END |