MINICERN/mathlib/gen/e/splas2.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:02:27  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10       SUBROUTINE SPLAS2(N,NC,NCX,NCY,NX,NY,MX,MY,KX,KY,NDIMC,NDIMZ,
  11      +                  XI,YI,ZI,KNOT,TX,TY,A,S,VT,W1,W2,LW,C,NERR)
  12
  13 #include "gen/imp64.inc"
  14       DIMENSION XI(*),YI(*),ZI(NDIMZ,*),TX(*),TY(*)
  15       DIMENSION A(N,*),S(*),VT(NC,*),W1(*),W2(*),C(NDIMC,*)
  16
  17 ************************************************************************
  18 *   NORBAS, VERSION: 15.03.1993
  19 ************************************************************************
  20 *
  21 *   THE SUBROUTINE  SPLAS2  IS USED BY  DSPAP2  FOR COMPUTING THE
  22 *   COEFFICIENTS
  23 *          C(I,J)   (I=1,...,NCX , J=1,...,NCY)
  24 *   OF A TWO-DIMENSIONAL POLYNOMIAL APPROXIMATION SPLINE  Z = S(X,Y)  IN
  25 *   REPRESENTATION OF NORMALIZED TWO-DIMENSIONAL B-SPLINES  B(I,J)(X,Y).
  26 *
  27 ************************************************************************
  28
  29       PARAMETER (Z0 = 0 , Z1 = 1 , Z2 = 2 , Z10 = 10 , HALF = Z1/Z2)
  30 *
  31 *   COMPUTE AN APPROXIMATION EPS0 TO THE RELATIVE MACHINE PRECISION
  32 *
  33       EPS0=Z1
  34    10 EPS0=EPS0/Z10
  35       IF(Z1+EPS0 .NE. Z1) GO TO 10
  36       EPS0=Z10*EPS0
  37 *
  38 *   COMPUTE KNOTS BY MEANS OF GIVEN DATA POINTS (IF KNOT = 1 OR 2)
  39 *
  40       IF (KNOT .EQ. 1) THEN
  41        CALL DSPKN2(KX,KY,MX,MY,XI(1),XI(NX),YI(1),YI(NY),TX,TY,NERR)
  42       ELSEIF(KNOT .EQ. 2) THEN
  43        DO 20 I=1,KX+1
  44        TX(I)=XI(1)
  45    20  TX(NCX+I)=XI(NX)
  46        DO 30 I=KX+2,NCX
  47    30  TX(I)=HALF*(XI(NX*(I-KX-2)/NCX+1)+XI(NX*I/NCX))
  48        DO 40 J=1,KY+1
  49        TY(J)=YI(1)
  50    40  TY(NCY+J)=YI(NY)
  51        DO 50 J=KY+2,NCY
  52    50  TY(J)=HALF*(YI(NY*(J-KY-2)/NCY+1)+YI(NY*J/NCY))
  53       ENDIF
  54 *
  55 *   COMPUTE MATRIX  A  AND SOLVE LINEAR LEAST SQUARES PROBLEM USING  SVD
  56 *
  57       DO 60  I=1,NX
  58       X=XI(I)
  59       DO 60  J=1,NY
  60       Y=YI(J)
  61       DO 60 IC=1,NCX
  62       DO 60 JC=1,NCY
  63       A((I-1)*NY+J,(IC-1)*NCY+JC)=
  64      +            DSPNB2(KX,KY,MX,MY,IC,JC,0,0,X,Y,TX,TY,NERR)
  65    60 CONTINUE
  66       CALL DGESVD('O','A',N,NC,A,N,S,W1,1,VT,NC,W2,LW,INFO)
  67       DO 70 I=1,NX
  68       DO 70 J=1,NY
  69    70 W1((I-1)*NY+J)=ZI(I,J)
  70       CALL DMMPY(NC,N,A(1,1),A(2,1),A(1,2),W1(1),W1(2),W2(1),W2(2))
  71       DO 80 J=1,NC
  72       IF(S(J) .GT. EPS0*S(1)) THEN
  73        W1(J)=W2(J)/S(J)
  74       ELSE
  75        W1(J)=Z0
  76       ENDIF
  77    80 CONTINUE
  78       CALL DMMPY(NC,NC,VT(1,1),VT(2,1),VT(1,2),W1(1),W1(2),W2(1),W2(2))
  79       DO 90 I=1,NCX
  80       DO 90 J=1,NCY
  81    90 C(I,J)=W2((I-1)*NCY+J)
  82       NERR=0
  83
  84       RETURN
  85       END
  86
  87
  88