MINICERN/mathlib/gen/e/dspap2.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 DSPAP2(KX,KY,MX,MY,NX,NY,XI,YI,ZI,NDIMZ,
  11      +                  KNOT,TX,TY,C,NDIMC,W,NW,NERR)
  12
  13 #include "gen/imp64.inc"
  14       DIMENSION XI(*),YI(*),ZI(NDIMZ,*),TX(*),TY(*),W(*),C(NDIMC,*)
  15       CHARACTER NAME*(*)
  16       CHARACTER*80 ERRTXT
  17       PARAMETER (NAME = 'DSPAP2')
  18
  19 ************************************************************************
  20 *   NORBAS, VERSION: 15.03.1993
  21 ************************************************************************
  22 *
  23 *   DSPAP2 COMPUTES THE COEFFICIENTS
  24 *          C(I,J)   (I=1,...,NCX , J=1,...,NCY)
  25 *   OF A TWO-DIMENSIONAL POLYNOMIAL APPROXIMATION SPLINE  Z = S(X,Y)  IN
  26 *   REPRESENTATION OF NORMALIZED TWO-DIMENSIONAL B-SPLINES  B(I,J)(X,Y)
  27 *
  28 *     S(X,Y) = SUMME(I=1,...,NCX)
  29 *              SUMME(J=1,...,NCY)  C(I,J) * B(I,J)(X,Y)
  30 *
  31 *   TO A USER SUPPLIED DATA SET
  32 *
  33 *     (XI(I),YI(J),ZI(I,J))       (I=1,...,NX , J=1,...,NY)
  34 *
  35 *   OF A FUNCTION  Z = F(X,Y) , I.E.
  36 *
  37 *     S(XI(I),YI(J)) = Z(I,J)     (I=1,...,NX , J=1,...,NY) .
  38 *
  39 *   THE TWO-DIMENSIONAL B-SPLINES  B(I,J)(X,Y)  ARE THE PRODUCT OF TWO
  40 *   ONE-DIMENSIONAL B-SPLINES  BX , BY
  41 *          B(I,J)(X,Y) = BX(I,KX)(X) * BY(J,KY)(Y)
  42 *   OF DEGREE  KX AND KY ( 0 <= KX , KY <= 25 )   WITH INDICES  I , J
  43 *   ( 1 <= I <= NX , 1 <= J <= NY ) OVER TWO SETS OF SPLINE-KNOTS
  44 *       TX(1),TX(2),...,TX(MX)     ( MX = NCX+KX+1 )
  45 *       TY(1),TY(2),...,TY(MY)     ( MY = NCY+KY+1 ) ,
  46 *   RESPECTIVELY.
  47 *   FOR FURTHER DETAILS TO THE ONE- AND TWO-DIMENSIONAL NORMALIZED
  48 *   B-SPLINES SEE THE COMMENTS TO  DSPNB1  AND  DSPNB2.
  49 *
  50 *   PARAMETERS:
  51 *
  52 *   NX    (INTEGER) NUMBER OF APPROXIMATION POINTS IN X-DIRECTION :
  53 *         XI(I) ,  I=1,...,NX .
  54 *   NY    (INTEGER) NUMBER OF APPROXIMATION POINTS IN Y-DIRECTION :
  55 *         YI(J) ,  J=1,...,NY .
  56 *   MX    (INTEGER) NUMBER OF KNOTS IN X-DIRECTION .
  57 *   MY    (INTEGER) NUMBER OF KNOTS IN Y-DIRECTION .
  58 *   KX    (INTEGER) DEGREE OF ONE-DIMENSIONAL B-SPLINES IN X-DIRECTION
  59 *         OVER THE SET OF KNOTS  TX.
  60 *   KY    (INTEGER) DEGREE OF ONE-DIMENSIONAL B-SPLINES IN Y-DIRECTION
  61 *         OVER THE SET OF KNOTS  TY.
  62 *   NDIMC (INTEGER) DECLARED FIRST DIMENSION OF ARRAY  C  IN THE
  63 *         CALLING PROGRAM, WITH  NDIMC >= (MX-KX-1) .
  64 *   NDIMZ (INTEGER) DECLARED FIRST DIMENSION OF ARRAY  ZI  IN THE
  65 *         CALLING PROGRAM, WITH  NDIMZ >= NX.
  66 *   XI    (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER NX .
  67 *         XI MUST CONTAIN THE APPROXIMATION POINTS IN X-DIRECTION IN
  68 *         ASCENDING ORDER, ON ENTRY.
  69 *   YI    (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER NY .
  70 *         YI MUST CONTAIN THE APPROXIMATION POINTS IN Y-DIRECTION IN
  71 *         ASCENDING ORDER, ON ENTRY.
  72 *   ZI    (DOUBLE PRECISION) ARRAY OF ORDER  (NDIMZ , >= NY) .
  73 *         ON ENTRY ZI MUST CONTAIN THE GIVEN FUNCTION VALUES
  74 *           Z(I,J)  AT THE APPROXIMATION POINTS  (X(I),Y(J))
  75 *                            ( I=1,...,NX , J=1,...,NY ).
  76 *   KNOT  (INTEGER) PARAMETER FOR STEERING THE CHOICE OF KNOTS.
  77 *         ON ENTRY:
  78 *         = 1 : KNOTS ARE COMPUTED BY  DSPAP2  IN THE FOLLOWING WAY:
  79 *               TX(J) = XI(1) ,                J = 1,...,KX+1
  80 *               TX(J) = XI(1)+(J-KX-1)*(XI(NX)-XI(1))/(NCX-KX) ,
  81 *                                              J = KX+2,...,NCX
  82 *               TX(NCX+J) = XI(NX) ,           J = 1,...,KX+1
  83 *         = 2 : KNOTS ARE COMPUTED BY  DSPAP2  IN THE FOLLOWING WAY:
  84 *               TX(J) = XI(1) ,                J = 1,...,KX+1
  85 *               TX(J) = (XI(NX*(J-KX-2)/NCX+1)+XI(NX*J/NCX))/2 ,
  86 *                                              J = KX+2,...,NCX
  87 *               TX(NCX+J) = XI(NX) ,           J = 1,...,KX+1
  88 *         OTHERWISE KNOTS ARE USER SUPPLIED. RECOMMENDED CHOICE :
  89 *               TX(1) <= ... <= TX(KX+1) <= XI(1)
  90 *               XI(1) < TX(KX+2) < ... < TX(NX) < XI(NX)
  91 *               XI(NX) <= TX(NX+1) <= ... <= TX(NX+KX+1)
  92 *         IN ALL CASES THE SAME CHOICE IS USED FOR KNOTS  TY  IN
  93 *         Y-DIRECTION .
  94 *   TX    (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER MX.
  95 *         IF THE INPUT VALUE OF THE PARAMETER KNOT IS  1 OR 2  THE
  96 *         KNOTS ARE COMPUTED BY  DSPAP2  AND THEY ARE GIVEN IN THE
  97 *         ARRAY TX, ON EXIT.
  98 *         IN THE OTHER CASES THE ARRAY  TX  MUST CONTAIN THE USER
  99 *         SUPPLIED KNOTS IN X-DIRECTION, ON ENTRY.
 100 *   TY    (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER MY.
 101 *         IF THE INPUT VALUE OF THE PARAMETER KNOT IS  1 OR 2  THE
 102 *         KNOTS ARE COMPUTED BY  DSPAP2  AND THEY ARE GIVEN IN THE
 103 *         ARRAY TY, ON EXIT.
 104 *         IN THE OTHER CASES THE ARRAY  TY  MUST CONTAIN THE USER
 105 *         SUPPLIED KNOTS IN Y-DIRECTION, ON ENTRY.
 106 *   W     (DOUBLE PRECISION) WORKING ARRAY OF AT LEAST ORDER NW .
 107 *   NW    (INTEGER) ORDER OF WORKING ARRAY W .
 108 *         NW >= N*(NC+6)+NC*(NC+1),
 109 *         WITH  N=NX*NY , NC=NCX*NCY=(MX-KX-1)*(MY-KY-1)
 110 *         FOR GOOD PERFORMANCE,  NW  SHOULD GENERALLY BE LARGER.
 111 *   C     (DOUBLE PRECISION) ARRAY OF ORDER (NDIMC, >= MY-KY-1).
 112 *         ON EXIT  C(I,J)  CONTAINS THE (I,J)-TH COEFFICIENT OF THE
 113 *         TWO-DIMENSIONAL B-SPLINE REPRESENTATION OF  S(X,Y) .
 114 *   NERR  (INTEGER) ERROR INDICATOR. ON EXIT:
 115 *         = 0: NO ERROR DETECTED
 116 *         = 1: AT LEAST ONE OF THE CONSTANTS  KX , KY , NX , NY , MX ,MY
 117 *              IS ILLEGAL
 118 *
 119 *   ERROR MESSAGES:
 120 *
 121 *   IF ONE OF THE FOLLOWING RELATION IS SATISFIED BY THE CHOSEN INPUT-
 122 *   PARAMETERS THE PROGRAM RETURNS, AND AN ERROR MESSAGE IS PRINTED:
 123 *     KX < 0         OR    KX > 25      OR
 124 *     KY < 0         OR    KY > 25      OR
 125 *     MX < 2*KX+2    OR    MY < 2*KY+2  OR
 126 *     NX < MX-KX-1   OR    NY < MY-KY-1 .
 127 *
 128 ************************************************************************
 129
 130       NERR=1
 131       IF(KX .LT. 0 .OR. KX .GT. 25) THEN
 132        WRITE(ERRTXT,101) 'KX',KX
 133        CALL MTLPRT(NAME,'E210.1',ERRTXT)
 134       ELSEIF(KY .LT. 0 .OR. KY .GT. 25) THEN
 135        WRITE(ERRTXT,101) 'KY',KY
 136        CALL MTLPRT(NAME,'E210.1',ERRTXT)
 137       ELSEIF(MX .LT. 2*KX+2) THEN
 138        WRITE(ERRTXT,101) 'MX',MX
 139        CALL MTLPRT(NAME,'E210.2',ERRTXT)
 140       ELSEIF(MY .LT. 2*KY+2) THEN
 141        WRITE(ERRTXT,101) 'MY',MY
 142        CALL MTLPRT(NAME,'E210.2',ERRTXT)
 143       ELSEIF(NX .LT. MX-KX-1) THEN
 144        WRITE(ERRTXT,101) 'NX',NX
 145        CALL MTLPRT(NAME,'E210.4',ERRTXT)
 146       ELSEIF(NY .LT. MY-KY-1) THEN
 147        WRITE(ERRTXT,101) 'NY',NY
 148        CALL MTLPRT(NAME,'E210.4',ERRTXT)
 149       ELSE
 150
 151        NCX=MX-KX-1
 152        NCY=MY-KY-1
 153        NC=NCX*NCY
 154        N=NX*NY
 155        M1=1
 156        M2=M1+N*NC
 157        M3=M2+NC
 158        M4=M3+NC*NC
 159        M5=M4+N
 160        LW=NW-M5+1
 161
 162        CALL SPLAS2(N,NC,NCX,NCY,NX,NY,MX,MY,KX,KY,NDIMC,NDIMZ,XI,YI,ZI,
 163      +             KNOT,TX,TY,W(M1),W(M2),W(M3),W(M4),W(M5),LW,C,NERR)
 164
 165       ENDIF
 166
 167       RETURN
 168
 169   101 FORMAT(1X,A5,' =',I6,'   NOT IN RANGE')
 170       END
 171
 172
 173