This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / e / cheb.F
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