This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / e / dspin2.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1996/04/01 15:02:25  mclareni
6 * Mathlib gen
7 *
8 *
9 #include "gen/pilot.h"
10       SUBROUTINE DSPIN2(KX,KY,NX,NY,XI,YI,ZI,NDIMZ,KNOT,TX,TY,C,NDIMC,
11      +                  W,IW,NERR)
12
13 #include "gen/imp64.inc"
14       DIMENSION XI(*),YI(*),ZI(NDIMZ,*),TX(*),TY(*)
15       DIMENSION W(*),IW(*),C(NDIMC,*)
16       CHARACTER NAME*(*)
17       CHARACTER*80 ERRTXT
18       PARAMETER (NAME = 'DSPIN2')
19
20
21 ************************************************************************
22 *   NORBAS, VERSION: 10.02.1993
23 ************************************************************************
24 *
25 *   DSPIN2 COMPUTES THE COEFFICIENTS
26 *          C(I,J)   (I=1,...,NX , J=1,...,NY)
27 *   OF A TWO-DIMENSIONAL POLYNOMIAL INTERPOLATION SPLINE  Z = S(X,Y)  IN
28 *   REPRESENTATION OF NORMALIZED TWO-DIMENSIONAL B-SPLINES  B(I,J)(X,Y)
29 *
30 *     S(X,Y) = SUMME(I=1,...,NX)
31 *              SUMME(J=1,...,NY)  C(I,J) * B(I,J)(X,Y)
32 *
33 *   TO A USER SUPPLIED DATA SET
34 *
35 *     (XI(I),YI(J),ZI(I,J))       (I=1,...,NX , J=1,...,NY)
36 *
37 *   OF A FUNCTION  Z = F(X,Y) , I.E.
38 *
39 *     S(XI(I),YI(J)) = Z(I,J)     (I=1,...,NX , J=1,...,NY) .
40 *
41 *   THE TWO-DIMENSIONAL B-SPLINES  B(I,J)(X,Y)  ARE THE PRODUCT OF TWO
42 *   ONE-DIMENSIONAL B-SPLINES  BX , BY
43 *          B(I,J)(X,Y) = BX(I,KX)(X) * BY(J,KY)(Y)
44 *   OF DEGREE  KX AND KY ( 0 <= KX , KY <= 25 )   WITH INDICES  I , J
45 *   ( 1 <= I <= NX , 1 <= J <= NY ) OVER TWO SETS OF SPLINE-KNOTS
46 *       TX(1),TX(2),...,TX(MX)     ( MX = NX+KX+1 )
47 *       TY(1),TY(2),...,TY(MY)     ( MY = NY+KY+1 ) ,
48 *   RESPECTIVELY.
49 *   FOR FURTHER DETAILS TO THE ONE- AND TWO-DIMENSIONAL NORMALIZED
50 *   B-SPLINES SEE THE COMMENTS TO  DSPNB1  AND  DSPNB2.
51 *
52 *   PARAMETERS:
53 *
54 *   NX    (INTEGER) NUMBER OF INTERPOLATION POINTS IN X-DIRECTION :
55 *         XI(I) ,  I=1,...,NX .
56 *   NY    (INTEGER) NUMBER OF INTERPOLATION POINTS IN Y-DIRECTION :
57 *         YI(J) ,  J=1,...,NY .
58 *   KX    (INTEGER) DEGREE OF ONE-DIMENSIONAL B-SPLINES IN X-DIRECTION
59 *         OVER THE SET OF KNOTS  TX,
60 *         WITH  KX <= NX-1  AND  0 <= KX <= 25 .
61 *   KY    (INTEGER) DEGREE OF ONE-DIMENSIONAL B-SPLINES IN Y-DIRECTION
62 *         OVER THE SET OF KNOTS  TY,
63 *         WITH  KY <= NY-1  AND  0 <= KY <= 25 .
64 *   NDIMC (INTEGER) DECLARED FIRST DIMENSION OF ARRAY  C  IN THE
65 *         CALLING PROGRAM, WITH  NDIMC >= NX .
66 *   NDIMZ (INTEGER) DECLARED FIRST DIMENSION OF ARRAY  ZI  IN THE
67 *         CALLING PROGRAM, WITH  NDIMZ >= NX .
68 *   XI    (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER NX .
69 *         XI MUST CONTAIN THE INTERPOLATION POINTS IN X-DIRECTION IN
70 *         ASCENDING ORDER, ON ENTRY.
71 *   YI    (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER NY .
72 *         YI MUST CONTAIN THE INTERPOLATION POINTS IN Y-DIRECTION IN
73 *         ASCENDING ORDER, ON ENTRY.
74 *   ZI    (DOUBLE PRECISION) ARRAY OF ORDER  (NDIMZ , >= NY) .
75 *         ON ENTRY ZI MUST CONTAIN THE GIVEN FUNCTION VALUES
76 *           Z(I,J)  AT THE INTERPOLATION POINTS  (X(I),Y(J))
77 *                            ( I=1,...,NX , J=1,...,NY ).
78 *   KNOT  (INTEGER) PARAMETER FOR STEERING THE CHOICE OF KNOTS.
79 *         ON ENTRY:
80 *         = 1 : KNOTS ARE COMPUTED BY  DSPIN2  IN THE FOLLOWING WAY:
81 *               TX(J) = XI(1) ,                J = 1,...,KX+1
82 *               TX(J) = XI(1)+(J-KX-1)*(XI(NX)-XI(1)) ,
83 *                                              J = KX+2,...,NX
84 *               TX(N+J) = XI(NX) ,             J = 1,...,KX+1
85 *         = 2 : KNOTS ARE COMPUTED BY  DSPIN2  IN THE FOLLOWING WAY:
86 *               TX(J) = XI(1) ,                J = 1,...,KX+1
87 *               TX(J) = (XI(J-KX-1)+XI(J))/2 , J = KX+2,...,NX
88 *               TX(N+J) = XI(NX) ,             J = 1,...,KX+1
89 *         OTHERWISE KNOTS ARE USER SUPPLIED. RECOMMENDED CHOICE :
90 *               T(1) <= ... <= T(K+1) <= XI(1)
91 *               XI(1) < T(K+2) < ... < T(N) < XI(N)
92 *               XI(N) <= T(N+1) <= ... <= T(N+K+1)
93 *         OTHERWISE KNOTS ARE USER SUPPLIED. RECOMMENDED CHOICE :
94 *               TX(1) <= ... <= TX(KX+1) <= XI(1)
95 *               XI(1) < TX(KX+2) < ... < TX(NX) < XI(NX)
96 *               XI(NX) <= TX(NX+1) <= ... <= TX(NX+KX+1)
97 *         IN ALL CASES THE SAME CHOICE IS USED FOR KNOTS  TY  IN
98 *         Y-DIRECTION .
99 *   TX    (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER MX.
100 *         IF THE INPUT VALUE OF THE PARAMETER KNOT IS  1 OR 2  THE
101 *         KNOTS ARE COMPUTED BY  DSPIN2  AND THEY ARE GIVEN IN THE
102 *         ARRAY TX, ON EXIT.
103 *         IN THE OTHER CASES THE ARRAY  TX  MUST CONTAIN THE USER
104 *         SUPPLIED KNOTS IN X-DIRECTION, ON ENTRY.
105 *   TY    (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER MY.
106 *         IF THE INPUT VALUE OF THE PARAMETER KNOT IS  1 OR 2  THE
107 *         KNOTS ARE COMPUTED BY  DSPIN2  AND THEY ARE GIVEN IN THE
108 *         ARRAY TY, ON EXIT.
109 *         IN THE OTHER CASES THE ARRAY  TY  MUST CONTAIN THE USER
110 *         SUPPLIED KNOTS IN Y-DIRECTION, ON ENTRY.
111 *   W     (DOUBLE PRECISION) WORKING ARRAY OF AT LEAST ORDER (L+1)*NINT,
112 *         WITH  NINT=NX*NY ,  K=KX*NY ,  L=3*K+1 .
113 *   IW    (INTEGER) WORKING ARRAY OF AT LEAST ORDER  NX*NY.
114 *   C     (DOUBLE PRECISION) ARRAY OF ORDER (NDIMC, >= NY).
115 *         ON EXIT  C(I,J)  CONTAINS THE (I,J)-TH COEFFICIENT OF THE
116 *         TWO-DIMENSIONAL B-SPLINE REPRESENTATION OF  S(X,Y) .
117 *   NERR  (INTEGER) ERROR INDICATOR. ON EXIT:
118 *         = 0: NO ERROR DETECTED
119 *         = 1: AT LEAST ONE OF THE CONSTANTS  KX , KY , NX , NY
120 *              IS ILLEGAL
121 *         = 2: THE LAPACK ROUTINES  DGBTRF , DGBTRS  COULD NOT SOLVE
122 *              THE LINEAR SYSTEM OF EQUATIONS WITH BAND-MATRIX FOR
123 *              COMPUTING C(1,1),...,C(NX,NY) . IT INDICATES THAT
124 *              A SOLUTION OF THE INTERPOLATION PROBLEM DOES NOT EXIST.
125 *              (ESPECIALLY, THE EXISTENCE OF A SOLUTION DEPENDS ON THE
126 *              SET OF KNOTS!)
127 *
128 *   ERROR MESSAGES:
129 *
130 *   IF ONE OF THE FOLLOWING RELATION IS SATISFIED BY THE CHOSEN INPUT-
131 *   PARAMETERS THE PROGRAM RETURNS, AND AN ERROR MESSAGE IS PRINTED:
132 *     KX < 0      OR    KX > 25    OR
133 *     KY < 0      OR    KY > 25    OR
134 *     NX < KX+1   OR    NY < KY+1 .
135 *
136 ************************************************************************
137
138
139       PARAMETER (Z1 = 1 , Z2 = 2 , HALF = Z1/Z2)
140
141       NERR=1
142       IF(KX .LT. 0 .OR. KX .GT. 25) THEN
143        WRITE(ERRTXT,101) 'KX',KX
144        CALL MTLPRT(NAME,'E210.1',ERRTXT)
145       ELSEIF(KY .LT. 0 .OR. KY .GT. 25) THEN
146        WRITE(ERRTXT,101) 'KY',KY
147        CALL MTLPRT(NAME,'E210.1',ERRTXT)
148       ELSEIF(NX .LT. KX+1) THEN
149        WRITE(ERRTXT,101) 'NX',NX
150        CALL MTLPRT(NAME,'E210.4',ERRTXT)
151       ELSEIF(NY .LT. KY+1) THEN
152        WRITE(ERRTXT,101) 'NY',NY
153        CALL MTLPRT(NAME,'E210.4',ERRTXT)
154       ELSE
155
156        K=KX*NY
157        L=3*K+1
158        NINT=NX*NY
159        MX=NX+KX+1
160        MY=NY+KY+1
161 *
162 *   COMPUTE KNOTS FROM INTERPOLATION POINTS ( IF KNOT=1 OR 2 )
163 *
164        IF (KNOT .EQ. 1) THEN
165         CALL DSPKN2(KX,KY,MX,MY,XI(1),XI(NX),YI(1),YI(NY),TX,TY,NERR)
166        ELSEIF (KNOT .EQ. 2) THEN
167         DO 10 I=1,KX+1
168         TX(I)=XI(1)
169    10   TX(NX+I)=XI(NX)
170         DO 20 I=KX+2,NX
171    20   TX(I)=HALF*(XI(I-KX-1)+XI(I))
172         DO 30 J=1,KY+1
173         TY(J)=YI(1)
174    30   TY(NY+J)=YI(NY)
175         DO 40 J=KY+2,NY
176    40   TY(J)=HALF*(YI(J-KY-1)+YI(J))
177        ENDIF
178 *
179 *   COMPUTE BAND-MATRIX  W
180 *
181 *   NUMBER OF SUB-, SUPER-DIAGONALS: KL=KU=(KX-1)*NY+KY-1<=KX*NY
182 *
183        DO 90 IC=1,NX
184        X=XI(IC)
185        DO 90 JC=1,NY
186        Y=YI(JC)
187        IJ=(IC-1)*NX+JC
188        DO 90 I=1,NX
189        DO 90 J=1,NY
190        JZ=(I-1)*NX+J
191        IZ=IJ+2*K+1-JZ
192        IF (K+1 .LE. IZ  .AND.  IZ .LE. 3*K+1)
193      +    W((JZ-1)*L+IZ)=DSPNB2(KX,KY,MX,MY,I,J,0,0,X,Y,TX,TY,NERR)
194    90  CONTINUE
195 *
196 *   SOLVE SYSTEM OF EQUATIONS FOR COMPUTING  C
197 *
198        DO 110 J=1,NINT
199   110  IW(J)=J
200        NERR=2
201        CALL DGBTRF(NINT,NINT,K,K,W,L,IW,INFO)
202        IF(INFO .NE. 0) RETURN
203        LW=L*NINT
204        DO 120 I=1,NX
205        DO 120 J=1,NY
206   120  W(LW+(I-1)*NY+J)=ZI(I,J)
207        CALL DGBTRS('N',NINT,K,K,1,W,L,IW,W(LW+1),NINT,INFO)
208        IF(INFO .NE. 0) RETURN
209        DO 130 I=1,NX
210        DO 130 J=1,NY
211   130  C(I,J)=W(LW+(I-1)*NY+J)
212        NERR=0
213       ENDIF
214
215       RETURN
216
217   101 FORMAT(1X,A5,' =',I6,'   NOT IN RANGE')
218       END
219
220
221