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

*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1996/04/01 15:02:25  mclareni
* Mathlib gen
*
*
#include "gen/pilot.h"
      SUBROUTINE DSPIN1(K,N,XI,YI,KNOT,T,C,W,IW,NERR)

#include "gen/imp64.inc"
      DIMENSION XI(*),YI(*),T(*),W(*),IW(*),C(*)
      CHARACTER NAME*(*)
      CHARACTER*80 ERRTXT
      PARAMETER (NAME = 'DSPIN1')

************************************************************************
*   NORBAS, VERSION: 10.02.1993
************************************************************************
*
*   DSPIN1 COMPUTES THE COEFFICIENTS  C(1),...,C(N)  OF A POLYNOMIAL
*   INTERPOLATION SPLINE  S(X)  IN B-SPLINE REPRESENTATION
*
*          S(X) =  SUMME(I=1,...,N)  C(I) * B(I,K)(X)
*
*   TO A USER SUPPLIED DATA SET
*
*          (XI(J),YI(J)) ,   J = 1,2,...,N  ,  N >= K+1
*
*   OF A FUNCTION  Y=F(X) , I.E.
*
*           S(XI(J)) = YI(J) , J = 1,2,...,N .
*
*   THE FUNCTIONS  B(I,K)(X)  ARE NORMALIZED B-SPLINES OF DEGREE  K
*   ( K <= N-1  AND  0<= K <= 25)  WITH INDEX  I  (1 <= I <= N) OVER A
*   SET OF SPLINE-KNOTS
*              T(1),T(2), ... ,T(M)    ( M = N+K+1 )
*   (KNOTS IN ASCENDING ORDER, WITH MULTIPLICITIES NOT GREATER
*   THAN  K+1).
*   FOR FURTHER DETAILS TO THE ONE-DIMENSIONAL NORMALIZED B-SPLINES SEE
*   THE COMMENTS TO  DSPNB1.
*
*   PARAMETERS:
*
*   N     (INTEGER) NUMBER OF INTERPOLATION POINTS .
*   K     (INTEGER) DEGREE (= ORDER - 1) OF B-SPLINES.
*   XI    (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER N .
*         XI MUST CONTAIN THE INTERPOLATION POINTS IN ASCENDING ORDER,
*         ON ENTRY.
*   YI    (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER N CONTAINING
*         THE FUNCTION VALUES YI(J), J=1,...,N,  ON ENTRY.
*   KNOT  (INTEGER) PARAMETER FOR STEERING THE CHOICE OF KNOTS.
*         ON ENTRY:
*         = 1 : KNOTS ARE COMPUTED BY  DSPIN1  IN THE FOLLOWING WAY:
*               T(J) = XI(1) ,                J = 1,...,K+1
*               T(J) = XI(1)+(J-K-1)*(XI(N)-XI(1)) ,
*                                             J = K+2,...,N
*               T(N+J) = XI(N) ,              J = 1,...,K+1
*         = 2 : KNOTS ARE COMPUTED BY  DSPIN1  IN THE FOLLOWING WAY:
*               T(J) = XI(1) ,                J = 1,...,K+1
*               T(J) = (XI(J-K-1)+XI(J))/2 ,  J = K+2,...,N
*               T(N+J) = XI(N) ,              J = 1,...,K+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)
*   T     (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER M .
*         IF THE INPUT VALUE OF THE PARAMETER KNOT IS  1 OR 2  THE
*         KNOTS ARE COMPUTED BY  DSPIN1  AND THEY ARE GIVEN IN THE
*         ARRAY T, ON EXIT.
*         IN THE OTHER CASES THE ARRAY  T  MUST CONTAIN THE USER
*         SUPPLIED KNOTS, ON ENTRY.
*   W     (DOUBLE PECISION) WORKING ARRAY OF AT LEAST ORDER
*         (3*K+1)*N .
*   IW    (INTEGER) WORKING ARRAY OF AT LEAST ORDER  N .
*   C     (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER  N . ON EXIT
*         C(1),...,C(N)  CONTAIN THE COEFFICIENTS OF THE B-SPLINE
*         REPRESENTATION OF  S(X).
*   NERR  (INTEGER) ERROR INDICATOR. ON EXIT:
*         = 0: NO ERROR DETECTED
*         = 1: AT LEAST ONE OF THE CONSTANTS  K , N  IS ILLEGAL
*         = 2: THE LAPACK ROUTINES  DGBTRF , DGBTRS  COULD NOT SOLVE
*              THE LINEAR SYSTEM OF EQUATIONS WITH BAND-MATRIX FOR
*              COMPUTING C(1),...,C(N) . 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:
*     K  < 0      OR    K > 25    OR
*     N < K+1 .
*
************************************************************************

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

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

       M=N+K+1
       L=3*K+1
*
*   COMPUTE KNOTS FROM INTERPOLATION POINTS (IF KNOT = 1 OR 2)
*
       IF (KNOT .EQ. 1) THEN
        CALL DSPKN1(K,M,XI(1),XI(N),T,NERR)
       ELSEIF (KNOT .EQ. 2) THEN
        DO 10 I=1,K+1
        T(I)=XI(1)
   10   T(N+I)=XI(N)
        DO 20 I=K+2,N
   20   T(I)=HALF*(XI(I-K-1)+XI(I))
       ENDIF
*
*   COMPUTE BAND-MATRIX  W
*
       DO 50 I=K+1,3*K+1
       DO 50 J=1,N
       IJ=I+J-2*K-1
       IF (1 .LE. IJ .AND. IJ .LE. N)
     +    W((J-1)*L+I)=DSPNB1(K,M,J,0,XI(IJ),T,NERR)
   50 CONTINUE
*
*   SOLVE SYSTEM OF EQUATIONS FOR COMPUTING  C
*
       DO 60 J=1,N
   60  IW(J)=J
       NERR=2
       CALL DGBTRF(N,N,K,K,W,L,IW,INFO)
       IF(INFO .NE. 0) RETURN
       CALL DVCPY(N,YI(1),YI(2),C(1),C(2))
       CALL DGBTRS('N',N,K,K,1,W,L,IW,C,N,INFO)
       IF(INFO .NE. 0) RETURN
       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 DSPIN1(K,N,XI,YI,KNOT,T,C,W,IW,NERR)
	11
	12	#include "gen/imp64.inc"
	13	DIMENSION XI(),YI(),T(),W(),IW(),C()
	14	CHARACTER NAME()
	15	CHARACTER*80 ERRTXT
	16	PARAMETER (NAME = 'DSPIN1')
	17
	18	************************************************************************
	19	* NORBAS, VERSION: 10.02.1993
	20	************************************************************************
	21	*
	22	* DSPIN1 COMPUTES THE COEFFICIENTS C(1),...,C(N) OF A POLYNOMIAL
	23	* INTERPOLATION SPLINE S(X) IN B-SPLINE REPRESENTATION
	24	*
	25	* S(X) = SUMME(I=1,...,N) C(I) * B(I,K)(X)
	26	*
	27	* TO A USER SUPPLIED DATA SET
	28	*
	29	* (XI(J),YI(J)) , J = 1,2,...,N , N >= K+1
	30	*
	31	* OF A FUNCTION Y=F(X) , I.E.
	32	*
	33	* S(XI(J)) = YI(J) , J = 1,2,...,N .
	34	*
	35	* THE FUNCTIONS B(I,K)(X) ARE NORMALIZED B-SPLINES OF DEGREE K
	36	* ( K <= N-1 AND 0<= K <= 25) WITH INDEX I (1 <= I <= N) OVER A
	37	* SET OF SPLINE-KNOTS
	38	* T(1),T(2), ... ,T(M) ( M = N+K+1 )
	39	* (KNOTS IN ASCENDING ORDER, WITH MULTIPLICITIES NOT GREATER
	40	* THAN K+1).
	41	* FOR FURTHER DETAILS TO THE ONE-DIMENSIONAL NORMALIZED B-SPLINES SEE
	42	* THE COMMENTS TO DSPNB1.
	43	*
	44	* PARAMETERS:
	45	*
	46	* N (INTEGER) NUMBER OF INTERPOLATION POINTS .
	47	* K (INTEGER) DEGREE (= ORDER - 1) OF B-SPLINES.
	48	* XI (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER N .
	49	* XI MUST CONTAIN THE INTERPOLATION POINTS IN ASCENDING ORDER,
	50	* ON ENTRY.
	51	* YI (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER N CONTAINING
	52	* THE FUNCTION VALUES YI(J), J=1,...,N, ON ENTRY.
	53	* KNOT (INTEGER) PARAMETER FOR STEERING THE CHOICE OF KNOTS.
	54	* ON ENTRY:
	55	* = 1 : KNOTS ARE COMPUTED BY DSPIN1 IN THE FOLLOWING WAY:
	56	* T(J) = XI(1) , J = 1,...,K+1
	57	* T(J) = XI(1)+(J-K-1)*(XI(N)-XI(1)) ,
	58	* J = K+2,...,N
	59	* T(N+J) = XI(N) , J = 1,...,K+1
	60	* = 2 : KNOTS ARE COMPUTED BY DSPIN1 IN THE FOLLOWING WAY:
	61	* T(J) = XI(1) , J = 1,...,K+1
	62	* T(J) = (XI(J-K-1)+XI(J))/2 , J = K+2,...,N
	63	* T(N+J) = XI(N) , J = 1,...,K+1
	64	* OTHERWISE KNOTS ARE USER SUPPLIED. RECOMMENDED CHOICE :
65	* T(1) <= ... <= T(K+1) <= XI(1)
66	* XI(1) < T(K+2) < ... < T(N) < XI(N)
67	* XI(N) <= T(N+1) <= ... <= T(N+K+1)
68	* T (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER M .
69	* IF THE INPUT VALUE OF THE PARAMETER KNOT IS 1 OR 2 THE
70	* KNOTS ARE COMPUTED BY DSPIN1 AND THEY ARE GIVEN IN THE
71	* ARRAY T, ON EXIT.
72	* IN THE OTHER CASES THE ARRAY T MUST CONTAIN THE USER
73	* SUPPLIED KNOTS, ON ENTRY.
74	* W (DOUBLE PECISION) WORKING ARRAY OF AT LEAST ORDER
75	* (3K+1)N .
76	* IW (INTEGER) WORKING ARRAY OF AT LEAST ORDER N .
77	* C (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER N . ON EXIT
78	* C(1),...,C(N) CONTAIN THE COEFFICIENTS OF THE B-SPLINE
79	* REPRESENTATION OF S(X).
80	* NERR (INTEGER) ERROR INDICATOR. ON EXIT:
81	* = 0: NO ERROR DETECTED
82	* = 1: AT LEAST ONE OF THE CONSTANTS K , N IS ILLEGAL
83	* = 2: THE LAPACK ROUTINES DGBTRF , DGBTRS COULD NOT SOLVE
84	* THE LINEAR SYSTEM OF EQUATIONS WITH BAND-MATRIX FOR
85	* COMPUTING C(1),...,C(N) . IT INDICATES THAT A SOLUTION
86	* OF THE INTERPOLATION PROBLEM DOES NOT EXIST.
87	* (ESPECIALLY, THE EXISTENCE OF A SOLUTION DEPENDS ON THE
88	* SET OF KNOTS!)
89	*
90	* ERROR MESSAGES:
91	*
92	* IF ONE OF THE FOLLOWING RELATION IS SATISFIED BY THE CHOSEN INPUT-
93	* PARAMETERS THE PROGRAM RETURNS, AND AN ERROR MESSAGE IS PRINTED:
94	* K < 0 OR K > 25 OR
95	* N < K+1 .
96	*
97	************************************************************************
98
99	PARAMETER (Z1 = 1 , Z2 = 2 , HALF = Z1/Z2)
100
101	NERR=1
102	IF(K .LT. 0 .OR. K .GT. 25) THEN
103	WRITE(ERRTXT,101) 'K',K
104	CALL MTLPRT(NAME,'E210.1',ERRTXT)
105	ELSEIF(N .LT. K+1) THEN
106	WRITE(ERRTXT,101) 'N',N
107	CALL MTLPRT(NAME,'E210.4',ERRTXT)
108	ELSE
109
110	M=N+K+1
111	L=3*K+1
112	*
113	* COMPUTE KNOTS FROM INTERPOLATION POINTS (IF KNOT = 1 OR 2)
114	*
115	IF (KNOT .EQ. 1) THEN
116	CALL DSPKN1(K,M,XI(1),XI(N),T,NERR)
117	ELSEIF (KNOT .EQ. 2) THEN
118	DO 10 I=1,K+1
119	T(I)=XI(1)
120	10 T(N+I)=XI(N)
121	DO 20 I=K+2,N
122	20 T(I)=HALF*(XI(I-K-1)+XI(I))
123	ENDIF
124	*
125	* COMPUTE BAND-MATRIX W
126	*
127	DO 50 I=K+1,3*K+1
128	DO 50 J=1,N
129	IJ=I+J-2*K-1
130	IF (1 .LE. IJ .AND. IJ .LE. N)
131	+ W((J-1)*L+I)=DSPNB1(K,M,J,0,XI(IJ),T,NERR)
132	50 CONTINUE
133	*
134	* SOLVE SYSTEM OF EQUATIONS FOR COMPUTING C
135	*
136	DO 60 J=1,N
137	60 IW(J)=J
138	NERR=2
139	CALL DGBTRF(N,N,K,K,W,L,IW,INFO)
140	IF(INFO .NE. 0) RETURN
141	CALL DVCPY(N,YI(1),YI(2),C(1),C(2))
142	CALL DGBTRS('N',N,K,K,1,W,L,IW,C,N,INFO)
143	IF(INFO .NE. 0) RETURN
144	NERR=0
145	ENDIF
146	RETURN
147
148	101 FORMAT(1X,A5,' =',I6,' NOT IN RANGE')
149	END
150
151
152