Mercurial > octave
view libcruft/villad/dif.f @ 8710:739141cde75a ss-3-1-52
fix typo in Array-f.cc
author | Jaroslav Hajek <highegg@gmail.com> |
---|---|
date | Mon, 09 Feb 2009 21:51:31 +0100 |
parents | 30c606bec7a8 |
children |
line wrap: on
line source
SUBROUTINE DIF ( NT, ROOT, DIF1, DIF2, DIF3 ) C INTEGER NT DOUBLE PRECISION ROOT(NT), DIF1(NT), DIF2(NT), DIF3(NT) C C*********************************************************************** C C SUBROUTINE DIF C C THIS ROUTINE IS NOT GIVEN SEPARATELY BY VILLADSEN AND MICHELSEN C BUT AS PART OF JCOBI C C DIF COMPUTES THE FIRST THREE DERIVATIVES OF THE NODE POLYNOMIAL C C N0 (ALPHA,BETA) N1 C P (X) = (X) * P (X) * (1 - X) C NT N C C AT THE INTERPOLATION POINTS. EACH OF THE PARAMETERS N0 AND N1 C MAY BE GIVEN THE VALUE 0 OR 1. NT = N + N0 + N1 C C THE VALUES OF ROOT MUST BE KNOWN BEFORE A CALL TO DIF IS POSSIBLE. C THEY MAY BE COMPUTED USING JCOBI. C C PARAMETER LIST: SEE THE SUBROUTINE JCOBI C C COMMON BLOCKS: NONE C C REQUIRED ROUTINES: VILERR C C*********************************************************************** C INTEGER I,J,IER DOUBLE PRECISION X,Y DOUBLE PRECISION ZERO,ONE,TWO,THREE LOGICAL LSTOP C PARAMETER ( ZERO = 0.0D+00, ONE = 1.0D+00, + TWO = 2.0D+00, THREE = 3.0D+00 ) C C -- ERROR CHECKING C IF (NT .LT. 1) THEN IER = 7 LSTOP = .TRUE. CALL VILERR(IER,LSTOP) ELSE END IF C C -- EVALUATE DERIVATIVES OF NODE POLYNOMIAL USING RECURSION FORMULAS C DO 40 I = 1,NT C X = ROOT(I) DIF1(I) = ONE DIF2(I) = ZERO DIF3(I) = ZERO C DO 30 J = 1,NT C IF (J .NE. I) THEN Y = X - ROOT(J) DIF3(I) = Y*DIF3(I) + THREE*DIF2(I) DIF2(I) = Y*DIF2(I) + TWO *DIF1(I) DIF1(I) = Y*DIF1(I) ELSE END IF C 30 CONTINUE 40 CONTINUE C RETURN END