[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / e / dspin2.F

*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1996/04/01 15:02:25  mclareni
* Mathlib gen
*
*
#include "gen/pilot.h"
      SUBROUTINE DSPIN2(KX,KY,NX,NY,XI,YI,ZI,NDIMZ,KNOT,TX,TY,C,NDIMC,
     +                  W,IW,NERR)

#include "gen/imp64.inc"
      DIMENSION XI(*),YI(*),ZI(NDIMZ,*),TX(*),TY(*)
      DIMENSION W(*),IW(*),C(NDIMC,*)
      CHARACTER NAME*(*)
      CHARACTER*80 ERRTXT
      PARAMETER (NAME = 'DSPIN2')


************************************************************************
*   NORBAS, VERSION: 10.02.1993
************************************************************************
*
*   DSPIN2 COMPUTES THE COEFFICIENTS
*          C(I,J)   (I=1,...,NX , J=1,...,NY)
*   OF A TWO-DIMENSIONAL POLYNOMIAL INTERPOLATION SPLINE  Z = S(X,Y)  IN
*   REPRESENTATION OF NORMALIZED TWO-DIMENSIONAL B-SPLINES  B(I,J)(X,Y)
*
*     S(X,Y) = SUMME(I=1,...,NX)
*              SUMME(J=1,...,NY)  C(I,J) * B(I,J)(X,Y)
*
*   TO A USER SUPPLIED DATA SET
*
*     (XI(I),YI(J),ZI(I,J))       (I=1,...,NX , J=1,...,NY)
*
*   OF A FUNCTION  Z = F(X,Y) , I.E.
*
*     S(XI(I),YI(J)) = Z(I,J)     (I=1,...,NX , J=1,...,NY) .
*
*   THE TWO-DIMENSIONAL B-SPLINES  B(I,J)(X,Y)  ARE THE PRODUCT OF TWO
*   ONE-DIMENSIONAL B-SPLINES  BX , BY
*          B(I,J)(X,Y) = BX(I,KX)(X) * BY(J,KY)(Y)
*   OF DEGREE  KX AND KY ( 0 <= KX , KY <= 25 )   WITH INDICES  I , J
*   ( 1 <= I <= NX , 1 <= J <= NY ) OVER TWO SETS OF SPLINE-KNOTS
*       TX(1),TX(2),...,TX(MX)     ( MX = NX+KX+1 )
*       TY(1),TY(2),...,TY(MY)     ( MY = NY+KY+1 ) ,
*   RESPECTIVELY.
*   FOR FURTHER DETAILS TO THE ONE- AND TWO-DIMENSIONAL NORMALIZED
*   B-SPLINES SEE THE COMMENTS TO  DSPNB1  AND  DSPNB2.
*
*   PARAMETERS:
*
*   NX    (INTEGER) NUMBER OF INTERPOLATION POINTS IN X-DIRECTION :
*         XI(I) ,  I=1,...,NX .
*   NY    (INTEGER) NUMBER OF INTERPOLATION POINTS IN Y-DIRECTION :
*         YI(J) ,  J=1,...,NY .
*   KX    (INTEGER) DEGREE OF ONE-DIMENSIONAL B-SPLINES IN X-DIRECTION
*         OVER THE SET OF KNOTS  TX,
*         WITH  KX <= NX-1  AND  0 <= KX <= 25 .
*   KY    (INTEGER) DEGREE OF ONE-DIMENSIONAL B-SPLINES IN Y-DIRECTION
*         OVER THE SET OF KNOTS  TY,
*         WITH  KY <= NY-1  AND  0 <= KY <= 25 .
*   NDIMC (INTEGER) DECLARED FIRST DIMENSION OF ARRAY  C  IN THE
*         CALLING PROGRAM, WITH  NDIMC >= NX .
*   NDIMZ (INTEGER) DECLARED FIRST DIMENSION OF ARRAY  ZI  IN THE
*         CALLING PROGRAM, WITH  NDIMZ >= NX .
*   XI    (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER NX .
*         XI MUST CONTAIN THE INTERPOLATION POINTS IN X-DIRECTION IN
*         ASCENDING ORDER, ON ENTRY.
*   YI    (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER NY .
*         YI MUST CONTAIN THE INTERPOLATION POINTS IN Y-DIRECTION IN
*         ASCENDING ORDER, ON ENTRY.
*   ZI    (DOUBLE PRECISION) ARRAY OF ORDER  (NDIMZ , >= NY) .
*         ON ENTRY ZI MUST CONTAIN THE GIVEN FUNCTION VALUES
*           Z(I,J)  AT THE INTERPOLATION POINTS  (X(I),Y(J))
*                            ( I=1,...,NX , J=1,...,NY ).
*   KNOT  (INTEGER) PARAMETER FOR STEERING THE CHOICE OF KNOTS.
*         ON ENTRY:
*         = 1 : KNOTS ARE COMPUTED BY  DSPIN2  IN THE FOLLOWING WAY:
*               TX(J) = XI(1) ,                J = 1,...,KX+1
*               TX(J) = XI(1)+(J-KX-1)*(XI(NX)-XI(1)) ,
*                                              J = KX+2,...,NX
*               TX(N+J) = XI(NX) ,             J = 1,...,KX+1
*         = 2 : KNOTS ARE COMPUTED BY  DSPIN2  IN THE FOLLOWING WAY:
*               TX(J) = XI(1) ,                J = 1,...,KX+1
*               TX(J) = (XI(J-KX-1)+XI(J))/2 , J = KX+2,...,NX
*               TX(N+J) = XI(NX) ,             J = 1,...,KX+1
*         OTHERWISE KNOTS ARE USER SUPPLIED. RECOMMENDED CHOICE :
*               T(1) <= ... <= T(K+1) <= XI(1)
*               XI(1) < T(K+2) < ... < T(N) < XI(N)
*               XI(N) <= T(N+1) <= ... <= T(N+K+1)
*         OTHERWISE KNOTS ARE USER SUPPLIED. RECOMMENDED CHOICE :
*               TX(1) <= ... <= TX(KX+1) <= XI(1)
*               XI(1) < TX(KX+2) < ... < TX(NX) < XI(NX)
*               XI(NX) <= TX(NX+1) <= ... <= TX(NX+KX+1)
*         IN ALL CASES THE SAME CHOICE IS USED FOR KNOTS  TY  IN
*         Y-DIRECTION .
*   TX    (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER MX.
*         IF THE INPUT VALUE OF THE PARAMETER KNOT IS  1 OR 2  THE
*         KNOTS ARE COMPUTED BY  DSPIN2  AND THEY ARE GIVEN IN THE
*         ARRAY TX, ON EXIT.
*         IN THE OTHER CASES THE ARRAY  TX  MUST CONTAIN THE USER
*         SUPPLIED KNOTS IN X-DIRECTION, ON ENTRY.
*   TY    (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER MY.
*         IF THE INPUT VALUE OF THE PARAMETER KNOT IS  1 OR 2  THE
*         KNOTS ARE COMPUTED BY  DSPIN2  AND THEY ARE GIVEN IN THE
*         ARRAY TY, ON EXIT.
*         IN THE OTHER CASES THE ARRAY  TY  MUST CONTAIN THE USER
*         SUPPLIED KNOTS IN Y-DIRECTION, ON ENTRY.
*   W     (DOUBLE PRECISION) WORKING ARRAY OF AT LEAST ORDER (L+1)*NINT,
*         WITH  NINT=NX*NY ,  K=KX*NY ,  L=3*K+1 .
*   IW    (INTEGER) WORKING ARRAY OF AT LEAST ORDER  NX*NY.
*   C     (DOUBLE PRECISION) ARRAY OF ORDER (NDIMC, >= NY).
*         ON EXIT  C(I,J)  CONTAINS THE (I,J)-TH COEFFICIENT OF THE
*         TWO-DIMENSIONAL B-SPLINE REPRESENTATION OF  S(X,Y) .
*   NERR  (INTEGER) ERROR INDICATOR. ON EXIT:
*         = 0: NO ERROR DETECTED
*         = 1: AT LEAST ONE OF THE CONSTANTS  KX , KY , NX , NY
*              IS ILLEGAL
*         = 2: THE LAPACK ROUTINES  DGBTRF , DGBTRS  COULD NOT SOLVE
*              THE LINEAR SYSTEM OF EQUATIONS WITH BAND-MATRIX FOR
*              COMPUTING C(1,1),...,C(NX,NY) . IT INDICATES THAT
*              A SOLUTION OF THE INTERPOLATION PROBLEM DOES NOT EXIST.
*              (ESPECIALLY, THE EXISTENCE OF A SOLUTION DEPENDS ON THE
*              SET OF KNOTS!)
*
*   ERROR MESSAGES:
*
*   IF ONE OF THE FOLLOWING RELATION IS SATISFIED BY THE CHOSEN INPUT-
*   PARAMETERS THE PROGRAM RETURNS, AND AN ERROR MESSAGE IS PRINTED:
*     KX < 0      OR    KX > 25    OR
*     KY < 0      OR    KY > 25    OR
*     NX < KX+1   OR    NY < KY+1 .
*
************************************************************************


      PARAMETER (Z1 = 1 , Z2 = 2 , HALF = Z1/Z2)

      NERR=1
      IF(KX .LT. 0 .OR. KX .GT. 25) THEN
       WRITE(ERRTXT,101) 'KX',KX
       CALL MTLPRT(NAME,'E210.1',ERRTXT)
      ELSEIF(KY .LT. 0 .OR. KY .GT. 25) THEN
       WRITE(ERRTXT,101) 'KY',KY
       CALL MTLPRT(NAME,'E210.1',ERRTXT)
      ELSEIF(NX .LT. KX+1) THEN
       WRITE(ERRTXT,101) 'NX',NX
       CALL MTLPRT(NAME,'E210.4',ERRTXT)
      ELSEIF(NY .LT. KY+1) THEN
       WRITE(ERRTXT,101) 'NY',NY
       CALL MTLPRT(NAME,'E210.4',ERRTXT)
      ELSE

       K=KX*NY
       L=3*K+1
       NINT=NX*NY
       MX=NX+KX+1
       MY=NY+KY+1
*
*   COMPUTE KNOTS FROM INTERPOLATION POINTS ( IF KNOT=1 OR 2 )
*
       IF (KNOT .EQ. 1) THEN
        CALL DSPKN2(KX,KY,MX,MY,XI(1),XI(NX),YI(1),YI(NY),TX,TY,NERR)
       ELSEIF (KNOT .EQ. 2) THEN
        DO 10 I=1,KX+1
        TX(I)=XI(1)
   10   TX(NX+I)=XI(NX)
        DO 20 I=KX+2,NX
   20   TX(I)=HALF*(XI(I-KX-1)+XI(I))
        DO 30 J=1,KY+1
        TY(J)=YI(1)
   30   TY(NY+J)=YI(NY)
        DO 40 J=KY+2,NY
   40   TY(J)=HALF*(YI(J-KY-1)+YI(J))
       ENDIF
*
*   COMPUTE BAND-MATRIX  W
*
*   NUMBER OF SUB-, SUPER-DIAGONALS: KL=KU=(KX-1)*NY+KY-1<=KX*NY
*
       DO 90 IC=1,NX
       X=XI(IC)
       DO 90 JC=1,NY
       Y=YI(JC)
       IJ=(IC-1)*NX+JC
       DO 90 I=1,NX
       DO 90 J=1,NY
       JZ=(I-1)*NX+J
       IZ=IJ+2*K+1-JZ
       IF (K+1 .LE. IZ  .AND.  IZ .LE. 3*K+1)
     +    W((JZ-1)*L+IZ)=DSPNB2(KX,KY,MX,MY,I,J,0,0,X,Y,TX,TY,NERR)
   90  CONTINUE
*
*   SOLVE SYSTEM OF EQUATIONS FOR COMPUTING  C
*
       DO 110 J=1,NINT
  110  IW(J)=J
       NERR=2
       CALL DGBTRF(NINT,NINT,K,K,W,L,IW,INFO)
       IF(INFO .NE. 0) RETURN
       LW=L*NINT
       DO 120 I=1,NX
       DO 120 J=1,NY
  120  W(LW+(I-1)*NY+J)=ZI(I,J)
       CALL DGBTRS('N',NINT,K,K,1,W,L,IW,W(LW+1),NINT,INFO)
       IF(INFO .NE. 0) RETURN
       DO 130 I=1,NX
       DO 130 J=1,NY
  130  C(I,J)=W(LW+(I-1)*NY+J)
       NERR=0
      ENDIF

      RETURN

  101 FORMAT(1X,A5,' =',I6,'   NOT IN RANGE')
      END
Commit	Line	Data
fe4da5cc	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=NXNY , K=KXNY , 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<=KXNY
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