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