Mercurial > octave-nkf
changeset 10603:d909c4c14b63
convert villad functions to C++
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Tue, 04 May 2010 13:00:00 -0400 |
| parents | 38eae0c3a003 |
| children | d7ff75c977e2 |
| files | ChangeLog ROADMAP libcruft/ChangeLog libcruft/Makefile.am libcruft/villad/dfopr.f libcruft/villad/dif.f libcruft/villad/intrp.f libcruft/villad/jcobi.f libcruft/villad/module.mk libcruft/villad/radau.f libcruft/villad/vilerr.f liboctave/ChangeLog liboctave/CollocWt.cc liboctave/CollocWt.h |
| diffstat | 14 files changed, 387 insertions(+), 929 deletions(-) [+] |
line wrap: on
line diff
--- a/ChangeLog Mon May 03 12:08:57 2010 -0700 +++ b/ChangeLog Tue May 04 13:00:00 2010 -0400 @@ -1,3 +1,7 @@ +2010-05-04 John W. Eaton <jwe@octave.org> + + * ROADMAP: Delete entry for villad. + 2010-04-26 Rik <octave@nomad.inbox5.com> * configure.ac: fix bug with shell HERE document introduced in
--- a/ROADMAP Mon May 03 12:08:57 2010 -0700 +++ b/ROADMAP Tue May 04 13:00:00 2010 -0400 @@ -27,7 +27,6 @@ ranlib * random number generators slatec-err * slatec error handling library slatec-fn * various special function subroutines - villad * subroutines for orthogonal collocation weights liboctave -- the C++ interfaces to the numerical libraries and various OS facilities
--- a/libcruft/ChangeLog Mon May 03 12:08:57 2010 -0700 +++ b/libcruft/ChangeLog Tue May 04 13:00:00 2010 -0400 @@ -1,3 +1,9 @@ +2010-05-04 John W. Eaton <jwe@octave.org> + + * villad/dfopr.f, villad/dif.f, villad/intrp.f, villad/jcobi.f, + villad/radau.f, villad/vilerr.f, villad/module.mk: Delete. + * Makefile.am: Don't include villad/module.mk. + 2010-04-11 Jaroslav Hajek <highegg@gmail.com> * blas-xtra/cmatm3.f, blas-xtra/zmatm3.f,
--- a/libcruft/Makefile.am Mon May 03 12:08:57 2010 -0700 +++ b/libcruft/Makefile.am Tue May 04 13:00:00 2010 -0400 @@ -71,7 +71,6 @@ include ranlib/module.mk include slatec-err/module.mk include slatec-fn/module.mk -include villad/module.mk cruft.def: $(libcruft_la_SOURCES) mkf77def chmod a+rx mkf77def
--- a/libcruft/villad/dfopr.f Mon May 03 12:08:57 2010 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,176 +0,0 @@ - SUBROUTINE DFOPR - + ( - + ND, N, N0, N1, I, ID, DIF1, DIF2, DIF3, ROOT, VECT - + ) - INTEGER ND, N, N0, N1, I, ID - DOUBLE PRECISION DIF1(ND), DIF2(ND), DIF3(ND), ROOT(ND), VECT(ND) -C -C*********************************************************************** -C -C VILLADSEN AND MICHELSEN, PAGES 133-134, 419 -C -C INPUT PARAMETERS: -C -C ND : THE DIMENSION OF THE VECTORS DIF1, DIF2, DIF3, AND ROOT -C -C N : THE DEGREE OF THE JACOBI POLYNOMIAL, (i.e. THE NUMBER -C OF INTERIOR INTERPOLATION POINTS) -C -C N0 : DETERMINES WHETHER X = 0 IS INCLUDED AS AN -C INTERPOLATION POINT -C -C N0 = 0 ==> X = 0 IS NOT INCLUDED -C N0 = 1 ==> X = 0 IS INCLUDED -C -C N1 : DETERMINES WHETHER X = 1 IS INCLUDED AS AN -C INTERPOLATION POINT -C -C N1 = 0 ==> X = 1 IS NOT INCLUDED -C N1 = 1 ==> X = 1 IS INCLUDED -C -C I : THE INDEX OF THE NODE FOR WHICH THE WEIGHTS ARE TO BE -C CALCULATED -C -C ID : INDICATOR -C -C ID = 1 ==> FIRST DERIVATIVE WEIGHTS ARE COMPUTED -C ID = 2 ==> SECOND DERIVATIVE WEIGHTS ARE COMPUTED -C ID = 3 ==> GAUSSIAN WEIGHTS ARE COMPUTED (IN THIS -C CASE, THE VALUE OF I IS IRRELEVANT) -C -C OUTPUT PARAMETERS: -C -C DIF1 : ONE DIMENSIONAL VECTOR CONTAINING THE FIRST DERIVATIVE -C OF THE NODE POLYNOMIAL AT THE ZEROS -C -C DIF2 : ONE DIMENSIONAL VECTOR CONTAINING THE SECOND DERIVATIVE -C OF THE NODE POLYNOMIAL AT THE ZEROS -C -C DIF3 : ONE DIMENSIONAL VECTOR CONTAINING THE THIRD DERIVATIVE -C OF THE NODE POLYNOMIAL AT THE ZEROS -C -C VECT : ONE DIMENSIONAL VECTOR OF COMPUTED WEIGHTS -C -C COMMON BLOCKS: NONE -C -C REQUIRED ROUTINES: VILERR -C -C*********************************************************************** -C - INTEGER J,NT,IER - DOUBLE PRECISION AX,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 ((N0 .NE. 0) .AND. (N0 .NE. 1)) THEN - IER = 1 - LSTOP = .TRUE. - CALL VILERR(IER,LSTOP) - ELSE - END IF -C - IF ((N1 .NE. 0) .AND. (N1 .NE. 1)) THEN - IER = 2 - LSTOP = .TRUE. - CALL VILERR(IER,LSTOP) - ELSE - END IF -C - IF (ND .LT. (N + N0 + N1)) THEN - IER = 3 - LSTOP = .TRUE. - CALL VILERR(IER,LSTOP) - ELSE - END IF -C - IF ((ID .NE. 1) .AND. (ID.NE. 2) .AND. (ID .NE. 3)) THEN - IER = 6 - LSTOP = .TRUE. - CALL VILERR(IER,LSTOP) - ELSE - END IF -C - IF (ID .NE. 3) THEN - IF (I .LT. 1) THEN - IER = 4 - LSTOP = .TRUE. - CALL VILERR(IER,LSTOP) - ELSE - END IF -C - IF (I .GT. (N + N0 + N1)) THEN - IER = 5 - LSTOP = .TRUE. - CALL VILERR(IER,LSTOP) - ELSE - END IF - ELSE - END IF -C - IF ((N + N0 + N1) .LT. 1) THEN - IER = 7 - LSTOP = .TRUE. - CALL VILERR(IER,LSTOP) - ELSE - END IF -C -C -- EVALUATE DISCRETIZATION MATRICES AND GAUSSIAN QUADRATURE -C -- WEIGHTS. QUADRATURE WEIGHTS ARE NORMALIZED TO SUM TO ONE. -C - NT = N + N0 + N1 -C - IF (ID .NE. 3) THEN - DO 20 J = 1,NT -C - IF (J .EQ. I) THEN - IF (ID .EQ. 1) THEN - VECT(I) = DIF2(I)/DIF1(I)/TWO - ELSE - VECT(I) = DIF3(I)/DIF1(I)/THREE - END IF - ELSE - Y = ROOT(I) - ROOT(J) - VECT(J) = DIF1(I)/DIF1(J)/Y - IF (ID .EQ. 2) THEN - VECT(J) = VECT(J)*(DIF2(I)/DIF1(I) - TWO/Y) - ELSE - END IF - END IF -C - 20 CONTINUE - ELSE - Y = ZERO -C - DO 25 J = 1,NT -C - X = ROOT(J) - AX = X*(ONE - X) -C - IF(N0 .EQ. 0) THEN - AX = AX/X/X - ELSE - END IF -C - IF(N1 .EQ. 0) THEN - AX = AX/(ONE - X)/(ONE - X) - ELSE - END IF -C - VECT(J) = AX/DIF1(J)**2 - Y = Y + VECT(J) -C - 25 CONTINUE -C - DO 60 J = 1,NT - VECT(J) = VECT(J)/Y - 60 CONTINUE -C - END IF -C - RETURN - END
--- a/libcruft/villad/dif.f Mon May 03 12:08:57 2010 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,74 +0,0 @@ - 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
--- a/libcruft/villad/intrp.f Mon May 03 12:08:57 2010 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,93 +0,0 @@ - SUBROUTINE INTRP ( ND, NT, X, ROOT, DIF1, XINTP ) -C - INTEGER ND, NT - DOUBLE PRECISION ROOT(ND), DIF1(ND), XINTP(ND) -C -C*********************************************************************** -C -C LAGRANGE INTERPOLATION -C -C VILLADSEN AND MICHELSEN, PAGES 132-133, 420 -C -C INPUT PARAMETERS: -C -C NT : THE TOTAL NUMBER OF INTERPOLATION POINTS FOR WHICH THE -C VALUE OF THE DEPENDENT VARIABLE Y IS KNOWN. NOTE: -C -C NT = N + N0 + N1 -C -C X : THE ABCISSA X WHERE Y(X) IS DESIRED -C -C ROOT : ONE DIMENSIONAL VECTOR CONTAINING ON EXIT THE -C N + N0 + N1 ZEROS OF THE NODE POLYNOMIAL USED IN THE -C INTERPOLATION ROUTINE -C -C DIF1 : ONE DIMENSIONAL VECTOR CONTAINING THE FIRST DERIVATIVE -C OF THE NODE POLYNOMIAL AT THE ZEROS -C -C OUTPUT PARAMETERS: -C -C XINTP : THE VECTOR OF INTERPOLATION WEIGHTS -C -C Y(X) IS GIVEN BY: -C -C NT -C Y(X) = SUM XINTRP(I) * Y(I) -C I=1 -C -C COMMON BLOCKS: NONE -C -C REQUIRED ROUTINES: VILERR -C -C*********************************************************************** -C - INTEGER I,IER - DOUBLE PRECISION POL,Y,X - DOUBLE PRECISION ZERO,ONE - LOGICAL LSTOP -C - PARAMETER ( ZERO = 0.0D+00, ONE = 1.0D+00 ) -C -C -- ERROR CHECKING -C - IF (ND .LT. NT) THEN - IER = 3 - LSTOP = .TRUE. - CALL VILERR(IER,LSTOP) - ELSE - END IF -C - IF (NT .LT. 1) THEN - IER = 7 - LSTOP = .TRUE. - CALL VILERR(IER,LSTOP) - ELSE - END IF -C -C -- EVALUATE LAGRANGIAN INTERPOLATION COEFFICIENTS -C - POL = ONE -C - DO 5 I = 1,NT -C - Y = X - ROOT(I) - XINTP(I) = ZERO -C - IF (Y .EQ. ZERO) THEN - XINTP(I) = ONE - ELSE - END IF -C - POL = POL*Y -C - 5 CONTINUE -C - IF (POL .NE. ZERO) THEN - DO 6 I = 1,NT - XINTP(I) = POL/DIF1(I)/(X - ROOT(I)) - 6 CONTINUE - ELSE - END IF -C - RETURN - END
--- a/libcruft/villad/jcobi.f Mon May 03 12:08:57 2010 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,240 +0,0 @@ -**************************************************************** -* -* The following routines (JCOBI, DIF, DFOPR, INTRP, AND RADAU) -* are the same as found in Villadsen, J. and M.L. Michelsen, -* Solution of Differential Equation Models by Polynomial -* Approximation, Prentice-Hall (1978) pages 418-420. -* -* Cosmetic changes (elimination of arithmetic IF statements, most -* GO TO statements, and indentation of program blocks) made by: -* -* John W. Eaton -* Department of Chemical Engineering -* The University of Texas at Austin -* Austin, Texas 78712 -* -* June 6, 1987 -* -* Some error checking additions also made on June 7, 1987 -* -* Further cosmetic changes made August 20, 1987 -* -************************************************************************ -* - SUBROUTINE JCOBI - + ( - + ND, N, N0, N1, ALPHA, BETA, DIF1, DIF2, DIF3, ROOT - + ) -C - INTEGER - + - + ND, N, N0, N1 -C - DOUBLE PRECISION - + - + ALPHA, BETA, DIF1(ND), DIF2(ND), DIF3(ND), ROOT(ND) -C -C*********************************************************************** -C -C VILLADSEN AND MICHELSEN, PAGES 131-132, 418 -C -C THIS SUBROUTINE COMPUTES THE ZEROS OF THE JACOBI POLYNOMIAL -C -C (ALPHA,BETA) -C P (X) -C N -C -C USE DIF (GIVEN BELOW) TO COMPUTE THE DERIVATIVES OF THE NODE -C POLYNOMIAL -C -C N0 (ALPHA,BETA) N1 -C P (X) = (X) * P (X) * (1 - X) -C NT N -C -C AT THE INTERPOLATION POINTS. -C -C INPUT PARAMETERS: -C -C ND : THE DIMENSION OF THE VECTORS DIF1, DIF2, DIF3, AND ROOT -C -C N : THE DEGREE OF THE JACOBI POLYNOMIAL, (i.e. THE NUMBER -C OF INTERIOR INTERPOLATION POINTS) -C -C N0 : DETERMINES WHETHER X = 0 IS INCLUDED AS AN -C INTERPOLATION POINT -C -C N0 = 0 ==> X = 0 IS NOT INCLUDED -C N0 = 1 ==> X = 0 IS INCLUDED -C -C N1 : DETERMINES WHETHER X = 1 IS INCLUDED AS AN -C INTERPOLATION POINT -C -C N1 = 0 ==> X = 1 IS NOT INCLUDED -C N1 = 1 ==> X = 1 IS INCLUDED -C -C ALPHA : THE VALUE OF ALPHA IN THE DESCRIPTION OF THE JACOBI -C POLYNOMIAL -C -C BETA : THE VALUE OF BETA IN THE DESCRIPTION OF THE JACOBI -C POLYNOMIAL -C -C FOR A MORE COMPLETE EXPLANATION OF ALPHA AN BETA, SEE VILLADSEN -C AND MICHELSEN, PAGES 57 TO 59 -C -C OUTPUT PARAMETERS: -C -C ROOT : ONE DIMENSIONAL VECTOR CONTAINING ON EXIT THE -C N + N0 + N1 ZEROS OF THE NODE POLYNOMIAL USED IN THE -C INTERPOLATION ROUTINE -C -C DIF1 : ONE DIMENSIONAL VECTOR CONTAINING THE FIRST DERIVATIVE -C OF THE NODE POLYNOMIAL AT THE ZEROS -C -C DIF2 : ONE DIMENSIONAL VECTOR CONTAINING THE SECOND DERIVATIVE -C OF THE NODE POLYNOMIAL AT THE ZEROS -C -C DIF3 : ONE DIMENSIONAL VECTOR CONTAINING THE THIRD DERIVATIVE -C OF THE NODE POLYNOMIAL AT THE ZEROS -C -C COMMON BLOCKS: NONE -C -C REQUIRED ROUTINES: VILERR, DIF -C -C*********************************************************************** -C - INTEGER I,J,NT,IER - DOUBLE PRECISION AB,AD,AP,Z1,Z,Y,X,XD,XN,XD1,XN1,XP,XP1,ZC - DOUBLE PRECISION ZERO,ONE,TWO - LOGICAL LSTOP -C - PARAMETER ( ZERO = 0.0D+00, ONE = 1.0D+00, TWO = 2.0D+00 ) -C -C -- ERROR CHECKING -C - IF ((N0 .NE. 0) .AND. (N0 .NE. 1)) THEN - IER = 1 - LSTOP = .TRUE. - CALL VILERR(IER,LSTOP) - ELSE - END IF -C - IF ((N1 .NE. 0) .AND. (N1 .NE. 1)) THEN - IER = 2 - LSTOP = .TRUE. - CALL VILERR(IER,LSTOP) - ELSE - END IF -C - IF (ND .LT. (N + N0 + N1)) THEN - IER = 3 - LSTOP = .TRUE. - CALL VILERR(IER,LSTOP) - ELSE - END IF -C - IF ((N + N0 + N1) .LT. 1) THEN - IER = 7 - LSTOP = .TRUE. - CALL VILERR(IER,LSTOP) - ELSE - END IF -C -C -- FIRST EVALUATION OF COEFFICIENTS IN RECURSION FORMULAS. -C -- RECURSION COEFFICIENTS ARE STORED IN DIF1 AND DIF2. -C - AB = ALPHA + BETA - AD = BETA - ALPHA - AP = BETA*ALPHA - DIF1(1) = (AD/(AB + TWO) + ONE)/TWO - DIF2(1) = ZERO -C - IF(N .GE. 2) THEN - DO 10 I = 2,N -C - Z1 = DBLE(I) - ONE - Z = AB + 2*Z1 - DIF1(I) = (AB*AD/Z/(Z + TWO) + ONE)/TWO -C - IF (I .EQ. 2) THEN - DIF2(I) = (AB + AP + Z1)/Z/Z/(Z + ONE) - ELSE - Z = Z*Z - Y = Z1*(AB + Z1) - Y = Y*(AP + Y) - DIF2(I) = Y/Z/(Z - ONE) - END IF -C - 10 CONTINUE - ELSE - END IF -C -C -- ROOT DETERMINATION BY NEWTON METHOD WITH SUPPRESSION OF -C -- PREVIOUSLY DETERMINED ROOTS -C - X = ZERO -C - DO 20 I = 1,N -C - 25 CONTINUE - XD = ZERO - XN = ONE - XD1 = ZERO - XN1 = ZERO -C - DO 30 J = 1,N - XP = (DIF1(J) - X)*XN - DIF2(J)*XD - XP1 = (DIF1(J) - X)*XN1 - DIF2(J)*XD1 - XN - XD = XN - XD1 = XN1 - XN = XP - XN1 = XP1 - 30 CONTINUE -C - ZC = ONE - Z = XN/XN1 -C - IF (I .NE. 1) THEN - DO 22 J = 2,I - ZC = ZC - Z/(X - ROOT(J-1)) - 22 CONTINUE - ELSE - END IF -C - Z = Z/ZC - X = X - Z -C - IF (DABS(Z) .GT. 1.D-09) THEN -C -C -- BACKWARD BRANCH -C - GO TO 25 - ELSE - END IF -C - ROOT(I) = X - X = X + 0.0001D0 -C - 20 CONTINUE -C -C -- ADD INTERPOLATION POINTS AT X = 0 AND/OR X = 1 -C - NT = N + N0 + N1 -C - IF (N0 .NE. 0) THEN - DO 31 I = 1,N - J = N + 1 - I - ROOT(J+1) = ROOT(J) - 31 CONTINUE - ROOT(1) = ZERO - ELSE - END IF -C - IF (N1 .EQ. 1) THEN - ROOT(NT) = ONE - ELSE - END IF -C - CALL DIF ( NT, ROOT, DIF1, DIF2, DIF3 ) -C - RETURN - END
--- a/libcruft/villad/module.mk Mon May 03 12:08:57 2010 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,9 +0,0 @@ -EXTRA_DIST += villad/module.mk - -libcruft_la_SOURCES += \ - villad/dfopr.f \ - villad/dif.f \ - villad/intrp.f \ - villad/jcobi.f \ - villad/radau.f \ - villad/vilerr.f
--- a/libcruft/villad/radau.f Mon May 03 12:08:57 2010 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,209 +0,0 @@ - SUBROUTINE RADAU - + ( - + ND, N, N0, N1, ID, ALPHA, BETA, ROOT, DIF1, VECT - + ) -C - INTEGER ND, N, N0, N1, ID - DOUBLE PRECISION ALPHA, BETA, ROOT(ND), DIF1(ND), VECT(ND) -C -C*********************************************************************** -C -C RADAU OR LOBATTO QUADRATURE -C -C VILLADSEN AND MICHELSEN, PAGES 133-135, 419 -C -C INPUT PARAMETERS: -C -C ND : THE DIMENSION OF THE VECTORS DIF1, DIF2, DIF3, AND ROOT -C -C N : THE DEGREE OF THE JACOBI POLYNOMIAL, (i.e. THE NUMBER -C OF INTERIOR INTERPOLATION POINTS) -C -C N0 : DETERMINES WHETHER X = 0 IS INCLUDED AS AN -C INTERPOLATION POINT -C -C N0 = 0 ==> X = 0 IS NOT INCLUDED -C N0 = 1 ==> X = 0 IS INCLUDED -C -C N1 : DETERMINES WHETHER X = 1 IS INCLUDED AS AN -C INTERPOLATION POINT -C -C N1 = 0 ==> X = 1 IS NOT INCLUDED -C N1 = 1 ==> X = 1 IS INCLUDED -C -C ID : INDICATOR -C -C ID = 1 ==> RADAU QUADRATURE WEIGHTS INCLUDING X = 1 -C ID = 2 ==> RADAU QUADRATURE WEIGHTS INCLUDING X = 0 -C ID = 3 ==> LOBATTO QUADRATURE WEIGHTS INCLUDING -C BOTH X = 0 AND X = 1 -C -C ALPHA : THE VALUE OF ALPHA IN THE DESCRIPTION OF THE JACOBI -C POLYNOMIAL -C -C BETA : THE VALUE OF BETA IN THE DESCRIPTION OF THE JACOBI -C POLYNOMIAL -C -C FOR A MORE COMPLETE EXPLANATION OF ALPHA AN BETA, SEE -C VILLADSEN AND MICHELSEN, PAGES 57 TO 59 -C -C ROOT : ONE DIMENSIONAL VECTOR CONTAINING ON EXIT THE -C N + N0 + N1 ZEROS OF THE NODE POLYNOMIAL USED IN THE -C INTERPOLATION ROUTINE -C -C DIF1 : ONE DIMENSIONAL VECTOR CONTAINING THE FIRST DERIVATIVE -C OF THE NODE POLYNOMIAL AT THE ZEROS -C -C THE NODE POLYNOMIAL IS GIVEN BY -C -C N0 (ALPHA',BETA') N1 -C P (X) = X * P (X) * (X - 1) -C NT N -C -C THE ARGUMENTS ALPHA' AND BETA' TO BE USED IN JCOBI FOR -C CALCULATION OF ROOT AND DIF1 DEPEND ON WHETHER X = 0 , X = 1 OR -C BOTH ARE USED AS EXTRA QUADRATURE POINTS. THUS: -C -C ID = 1: ALPHA' = ALPHA + 1, BETA' = BETA -C ID = 2: ALPHA' = ALPHA , BETA' = BETA + 1 -C ID = 3: ALPHA' = ALPHA + 1, BETA' = BETA + 1 -C -C NOTE: -C -C ID = 1 REQUIRES THAT N0 = 0 OR 1, N1 = 1 -C ID = 2 REQUIRES THAT N0 = 1 , N1 = 0 OR 1 -C ID = 3 REQUIRES THAT N0 = 1 , N1 = 1 -C -C OUTPUT PARAMETERS: -C -C VECT : VECTOR OF THE NT COMPUTED QUADRATURE WEIGHTS, -C NORMALIZED SUCH THAT -C -C NT -C SUM VECT(I) = 1 -C I=1 -C -C FOR A MORE COMPLETE EXPLANATION SEE VILLADSEN AND -C MICHELSEN, PAGES 133 TO 135 -C -C COMMON BLOCKS: NONE -C -C REQUIRED ROUTINES: VILERR -C -C*********************************************************************** -C - INTEGER I,NT,IER - DOUBLE PRECISION AX,S,X - DOUBLE PRECISION ZERO,ONE - LOGICAL LSTOP -C - PARAMETER ( ZERO = 0.0D+00, ONE = 1.0D+00 ) -C -C -- ERROR CHECKING -C - IF ((N0 .NE. 0) .AND. (N0 .NE. 1)) THEN - IER = 1 - LSTOP = .TRUE. - CALL VILERR(IER,LSTOP) - ELSE - END IF -C - IF ((N1 .NE. 0) .AND. (N1 .NE. 1)) THEN - IER = 2 - LSTOP = .TRUE. - CALL VILERR(IER,LSTOP) - ELSE - END IF -C - IF (ND .LT. (N + N0 + N1)) THEN - IER = 3 - LSTOP = .TRUE. - CALL VILERR(IER,LSTOP) - ELSE - END IF -C - IF ((N + N0 + N1) .LT. 1) THEN - IER = 7 - LSTOP = .TRUE. - CALL VILERR(IER,LSTOP) - ELSE - END IF -C - IF ((ID .NE. 1) .AND. (ID.NE. 2) .AND. (ID .NE. 3)) THEN - IER = 8 - LSTOP = .TRUE. - CALL VILERR(IER,LSTOP) - ELSE - END IF -C - IF ((ID .EQ. 1) .AND. (N1 .NE. 1)) THEN - IER = 9 - LSTOP = .TRUE. - CALL VILERR(IER,LSTOP) - ELSE - END IF -C - IF ((ID .EQ. 2) .AND. (N0 .NE. 1)) THEN - IER = 10 - LSTOP = .TRUE. - CALL VILERR(IER,LSTOP) - ELSE - END IF -C - IF ((ID .EQ. 3) .AND. ((N0 .NE. 1) .OR. (N1 .NE. 1))) THEN - IER = 11 - LSTOP = .TRUE. - CALL VILERR(IER,LSTOP) - ELSE - END IF -C -C -- EVALUATE RADAU OR LOBATTO QUADRATURE WEIGHTS -C - S = ZERO - NT = N + N0 + N1 -C - DO 40 I = 1,NT -C - X = ROOT(I) -C - IF (ID .EQ. 1) THEN - AX = X - IF (N0 .EQ. 0) THEN - AX = ONE/AX - ELSE - END IF - ELSE IF (ID .EQ. 2) THEN - AX = ONE - X - IF (N1 .EQ. 0) THEN - AX = ONE/AX - ELSE - END IF - ELSE IF (ID .EQ. 3) THEN - AX = ONE - ELSE - END IF -C - VECT(I) = AX/DIF1(I)**2 -C - 40 CONTINUE -C - IF (ID .NE. 2) THEN - VECT(NT) = VECT(NT)/(ONE + ALPHA) - ELSE - END IF -C - IF (ID .GT. 1) THEN - VECT(1) = VECT( 1)/(ONE + BETA) - ELSE - END IF -C - DO 50 I = 1,NT - S = S + VECT(I) - 50 CONTINUE -C - DO 60 I = 1,NT - VECT(I) = VECT(I)/S - 60 CONTINUE -C - RETURN - END
--- a/libcruft/villad/vilerr.f Mon May 03 12:08:57 2010 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,89 +0,0 @@ - SUBROUTINE VILERR ( IER, LSTOP ) -C - INTEGER IER - LOGICAL LSTOP -C -C*********************************************************************** -C -C THIS SUBROUTINE HANDLES ERRORS FOR THE SUBROUTINES JCOBI, DFOPR, -C INTRP, AND RADAU GIVEN BY VILLADSEN AND MICHELSEN. -C -C PARAMETER LIST: -C -C IER : ERROR NUMBER -C LSTOP : LOGICAL FLAG -C -C LSTOP = .TRUE. ==> FATAL ERROR, PROGRAM TERMINATION -C LSTOP = .FALSE. ==> WARNING ERROR, NORMAL RETURN -C -C COMMON BLOCKS: NONE -C -C REQUIRED ROUTINES: NONE -C -C*********************************************************************** -C -C -- BEGIN -C - IF ( IER .EQ. 1) THEN -C - WRITE(*,*) '** VILERR : Illegal value for N0 ' -C - ELSE IF ( IER .EQ. 2) THEN -C - WRITE(*,*) '** VILERR : Illegal value for N1 ' -C - ELSE IF ( IER .EQ. 3 ) THEN -C - WRITE(*,*) '** VILERR : Insufficient dimension for problem ' -C - ELSE IF ( IER .EQ. 4 ) THEN -C - WRITE(*,*) '** VILERR : Index less than zero in DFOPR ' -C - ELSE IF ( IER .EQ. 5 ) THEN -C - WRITE(*,*) '** VILERR : Index greater than NTOTAL in DFOPR ' -C - ELSE IF ( IER .EQ. 6 ) THEN -C - WRITE(*,*) '** VILERR : Illegal ID in DFOPR ' -C - ELSE IF ( IER .EQ. 7 ) THEN -C - WRITE(*,*) '** VILERR : Number of interpolation points ' - WRITE(*,*) ' less than 1 ' -C - ELSE IF ( IER .EQ. 8 ) THEN -C - WRITE(*,*) '** VILERR : Illegal ID in RADAU ' -C - ELSE IF ( IER .EQ. 9 ) THEN -C - WRITE(*,*) '** VILERR : ID = 1 but N1 not equal to 1 in RADAU ' -C - ELSE IF ( IER .EQ. 10 ) THEN -C - WRITE(*,*) '** VILERR : ID = 2 but N0 not equal to 1 in RADAU ' -C - ELSE IF ( IER .EQ. 11 ) THEN -C - WRITE(*,*) '** VILERR : ID = 3 but N0 not equal to 1 or ' - WRITE(*,*) ' N1 not equal to 1 in RADAU ' -C - ELSE -C - WRITE(*,*) 'UNRECOGNIZED ERROR FLAG SET FOR VILERR ' -C - END IF -C - IF ( LSTOP ) THEN -C -C -- PROGRAM EXECUTION TERMINATES HERE -C - CALL XSTOPX (' ') -C - ELSE - END IF -C - RETURN - END
--- a/liboctave/ChangeLog Mon May 03 12:08:57 2010 -0700 +++ b/liboctave/ChangeLog Tue May 04 13:00:00 2010 -0400 @@ -1,3 +1,13 @@ +2010-05-04 John W. Eaton <jwe@octave.org> + + * CollocWt.cc (diff, jcobi, dfopr): New functions, based on + Fortran functions in libcruft/villad. + (jcobi): Handle iteration failure at large N. + (CollocWt::init): Call them instead of Fortran code. + * CollocWt.h (CollocWt::initialized): Declare as bool, not int. + Change all uses. + Addresses bug #29473. + 2010-05-03 Jaroslav Hajek <highegg@gmail.com> * dbleSVD.h (SVD::driver): New enum.
--- a/liboctave/CollocWt.cc Mon May 03 12:08:57 2010 -0700 +++ b/liboctave/CollocWt.cc Tue May 04 13:00:00 2010 -0400 @@ -27,21 +27,348 @@ #include <iostream> +#include <cfloat> + #include "CollocWt.h" #include "f77-fcn.h" #include "lo-error.h" -extern "C" +// The following routines jcobi, dif, and dfopr are based on the code +// found in Villadsen, J. and M. L. Michelsen, Solution of Differential +// Equation Models by Polynomial Approximation, Prentice-Hall (1978) +// pages 418-420. +// +// Translated to C++ by jwe. + +// Compute the first three derivatives of the node polynomial. +// +// n0 (alpha,beta) n1 +// p (x) = (x) * p (x) * (1 - x) +// nt n +// +// at the interpolation points. Each of the parameters n0 and n1 +// may be given the value 0 or 1. The total number of points +// nt = n + n0 + n1 +// +// The values of root must be known before a call to dif is possible. +// They may be computed using jcobi. + +static void +dif (octave_idx_type nt, double *root, double *dif1, double *dif2, + double *dif3) +{ + // Evaluate derivatives of node polynomial using recursion formulas. + + for (octave_idx_type i = 0; i < nt; i++) + { + double x = root[i]; + + dif1[i] = 1.0; + dif2[i] = 0.0; + dif3[i] = 0.0; + + for (octave_idx_type j = 0; j < nt; j++) + { + if (j != i) + { + double y = x - root[j]; + + dif3[i] = y * dif3[i] + 3.0 * dif2[i]; + dif2[i] = y * dif2[i] + 2.0 * dif1[i]; + dif1[i] = y * dif1[i]; + } + } + } +} + +// Compute the zeros of the Jacobi polynomial. +// +// (alpha,beta) +// p (x) +// n +// +// Use dif to compute the derivatives of the node +// polynomial +// +// n0 (alpha,beta) n1 +// p (x) = (x) * p (x) * (1 - x) +// nt n +// +// at the interpolation points. +// +// See Villadsen and Michelsen, pages 131-132 and 418. +// +// Input parameters: +// +// nd : the dimension of the vectors dif1, dif2, dif3, and root +// +// n : the degree of the jacobi polynomial, (i.e. the number +// of interior interpolation points) +// +// n0 : determines whether x = 0 is included as an +// interpolation point +// +// n0 = 0 ==> x = 0 is not included +// n0 = 1 ==> x = 0 is included +// +// n1 : determines whether x = 1 is included as an +// interpolation point +// +// n1 = 0 ==> x = 1 is not included +// n1 = 1 ==> x = 1 is included +// +// alpha : the value of alpha in the description of the jacobi +// polynomial +// +// beta : the value of beta in the description of the jacobi +// polynomial +// +// For a more complete explanation of alpha an beta, see Villadsen +// and Michelsen, pages 57 to 59. +// +// Output parameters: +// +// root : one dimensional vector containing on exit the +// n + n0 + n1 zeros of the node polynomial used in the +// interpolation routine +// +// dif1 : one dimensional vector containing the first derivative +// of the node polynomial at the zeros +// +// dif2 : one dimensional vector containing the second derivative +// of the node polynomial at the zeros +// +// dif3 : one dimensional vector containing the third derivative +// of the node polynomial at the zeros + +static bool +jcobi (octave_idx_type n, octave_idx_type n0, octave_idx_type n1, + double alpha, double beta, double *dif1, double *dif2, + double *dif3, double *root) { - F77_RET_T - F77_FUNC (jcobi, JCOBI) (octave_idx_type&, octave_idx_type&, octave_idx_type&, octave_idx_type&, double&, - double&, double*, double*, double*, - double*); + assert (n0 == 0 || n0 == 1); + assert (n1 == 0 || n1 == 1); + + octave_idx_type nt = n + n0 + n1; + + assert (nt > 1); + +// -- first evaluation of coefficients in recursion formulas. +// -- recursion coefficients are stored in dif1 and dif2. + + double ab = alpha + beta; + double ad = beta - alpha; + double ap = beta * alpha; + + dif1[0] = (ad / (ab + 2.0) + 1.0) / 2.0; + dif2[0] = 0.0; + + if (n >= 2) + { + for (octave_idx_type i = 1; i < n; i++) + { + double z1 = i; + double z = ab + 2 * z1; + + dif1[i] = (ab * ad / z / (z + 2.0) + 1.0) / 2.0; + + if (i == 1) + dif2[i] = (ab + ap + z1) / z / z / (z + 1.0); + else + { + z = z * z; + double y = z1 * (ab + z1); + y = y * (ap + y); + dif2[i] = y / z / (z - 1.0); + } + } + } + + // Root determination by Newton method with suppression of previously + // determined roots. + + double x = 0.0; + + for (octave_idx_type i = 0; i < n; i++) + { + bool done = false; + + int k = 0; + + while (! done) + { + double xd = 0.0; + double xn = 1.0; + double xd1 = 0.0; + double xn1 = 0.0; + + for (octave_idx_type j = 0; j < n; j++) + { + double xp = (dif1[j] - x) * xn - dif2[j] * xd; + double xp1 = (dif1[j] - x) * xn1 - dif2[j] * xd1 - xn; + + xd = xn; + xd1 = xn1; + xn = xp; + xn1 = xp1; + } + + double zc = 1.0; + double z = xn / xn1; + + if (i != 0) + { + for (octave_idx_type j = 1; j <= i; j++) + zc = zc - z / (x - root[j-1]); + } + + z = z / zc; + x = x - z; + + // Famous last words: 100 iterations should be more than + // enough in all cases. + + if (++k > 100 || xisnan (z)) + return false; + + if (std::abs (z) <= 100 * DBL_EPSILON) + done = true; + } + + root[i] = x; + x = x + sqrt (DBL_EPSILON); + } + + // Add interpolation points at x = 0 and/or x = 1. + + if (n0 != 0) + { + for (octave_idx_type i = n; i > 0; i--) + root[i] = root[i-1]; + + root[0] = 0.0; + } + + if (n1 != 0) + root[nt-1] = 1.0; + + dif (nt, root, dif1, dif2, dif3); + + return true; +} - F77_RET_T - F77_FUNC (dfopr, DFOPR) (octave_idx_type&, octave_idx_type&, octave_idx_type&, octave_idx_type&, octave_idx_type&, octave_idx_type&, - double*, double*, double*, double*, - double*); +// Compute derivative weights for orthogonal collocation. +// +// See Villadsen and Michelsen, pages 133-134, 419. +// +// Input parameters: +// +// nd : the dimension of the vectors dif1, dif2, dif3, and root +// +// n : the degree of the jacobi polynomial, (i.e. the number +// of interior interpolation points) +// +// n0 : determines whether x = 0 is included as an +// interpolation point +// +// n0 = 0 ==> x = 0 is not included +// n0 = 1 ==> x = 0 is included +// +// n1 : determines whether x = 1 is included as an +// interpolation point +// +// n1 = 0 ==> x = 1 is not included +// n1 = 1 ==> x = 1 is included +// +// i : the index of the node for which the weights are to be +// calculated +// +// id : indicator +// +// id = 1 ==> first derivative weights are computed +// id = 2 ==> second derivative weights are computed +// id = 3 ==> gaussian weights are computed (in this +// case, the value of i is irrelevant) +// +// Output parameters: +// +// dif1 : one dimensional vector containing the first derivative +// of the node polynomial at the zeros +// +// dif2 : one dimensional vector containing the second derivative +// of the node polynomial at the zeros +// +// dif3 : one dimensional vector containing the third derivative +// of the node polynomial at the zeros +// +// vect : one dimensional vector of computed weights + +static void +dfopr (octave_idx_type n, octave_idx_type n0, octave_idx_type n1, + octave_idx_type i, octave_idx_type id, double *dif1, + double *dif2, double *dif3, double *root, double *vect) +{ + assert (n0 == 0 || n0 == 1); + assert (n1 == 0 || n1 == 1); + + octave_idx_type nt = n + n0 + n1; + + assert (nt > 1); + + assert (id == 1 || id == 2 || id == 3); + + if (id != 3) + assert (i >= 0 && i < nt); + + // Evaluate discretization matrices and Gaussian quadrature weights. + // Quadrature weights are normalized to sum to one. + + if (id != 3) + { + for (octave_idx_type j = 0; j < nt; j++) + { + if (j == i) + { + if (id == 1) + vect[i] = dif2[i] / dif1[i] / 2.0; + else + vect[i] = dif3[i] / dif1[i] / 3.0; + } + else + { + double y = root[i] - root[j]; + + vect[j] = dif1[i] / dif1[j] / y; + + if (id == 2) + vect[j] = vect[j] * (dif2[i] / dif1[i] - 2.0 / y); + } + } + } + else + { + double y = 0.0; + + for (octave_idx_type j = 0; j < nt; j++) + { + double x = root[j]; + + double ax = x * (1.0 - x); + + if (n0 == 0) + ax = ax / x / x; + + if (n1 == 0) + ax = ax / (1.0 - x) / (1.0 - x); + + vect[j] = ax / (dif1[j] * dif1[j]); + + y = y + vect[j]; + } + + for (octave_idx_type j = 0; j < nt; j++) + vect[j] = vect[j] / y; + } } // Error handling. @@ -123,41 +450,41 @@ // Compute roots. - F77_FUNC (jcobi, JCOBI) (nt, n, inc_left, inc_right, Alpha, Beta, - pdif1, pdif2, pdif3, pr); + if (! jcobi (n, inc_left, inc_right, Alpha, Beta, pdif1, pdif2, pdif3, pr)) + { + error ("jcobi: newton iteration failed"); + return; + } octave_idx_type id; // First derivative weights. id = 1; - for (octave_idx_type i = 1; i <= nt; i++) + for (octave_idx_type i = 0; i < nt; i++) { - F77_FUNC (dfopr, DFOPR) (nt, n, inc_left, inc_right, i, id, pdif1, - pdif2, pdif3, pr, pvect); + dfopr (n, inc_left, inc_right, i, id, pdif1, pdif2, pdif3, pr, pvect); for (octave_idx_type j = 0; j < nt; j++) - A (i-1, j) = vect.elem (j); + A(i,j) = vect(j); } // Second derivative weights. id = 2; - for (octave_idx_type i = 1; i <= nt; i++) + for (octave_idx_type i = 0; i < nt; i++) { - F77_FUNC (dfopr, DFOPR) (nt, n, inc_left, inc_right, i, id, pdif1, - pdif2, pdif3, pr, pvect); + dfopr (n, inc_left, inc_right, i, id, pdif1, pdif2, pdif3, pr, pvect); for (octave_idx_type j = 0; j < nt; j++) - B (i-1, j) = vect.elem (j); + B(i,j) = vect(j); } // Gaussian quadrature weights. id = 3; double *pq = q.fortran_vec (); - F77_FUNC (dfopr, DFOPR) (nt, n, inc_left, inc_right, id, id, pdif1, - pdif2, pdif3, pr, pq); + dfopr (n, inc_left, inc_right, id, id, pdif1, pdif2, pdif3, pr, pq); initialized = 1; }
--- a/liboctave/CollocWt.h Mon May 03 12:08:57 2010 -0700 +++ b/liboctave/CollocWt.h Tue May 04 13:00:00 2010 -0400 @@ -37,24 +37,27 @@ CollocWt (void) : n (0), inc_left (0), inc_right (0), lb (0.0), rb (1.0), - Alpha (0.0), Beta (0.0), r (), q (), A (), B (), initialized (0) { } + Alpha (0.0), Beta (0.0), r (), q (), A (), B (), initialized (false) { } CollocWt (octave_idx_type nc, octave_idx_type il, octave_idx_type ir) : n (nc), inc_left (il), inc_right (ir), lb (0.0), rb (1.0), - Alpha (0.0), Beta (0.0), r (), q (), A (), B (), initialized (0) { } + Alpha (0.0), Beta (0.0), r (), q (), A (), B (), initialized (false) { } - CollocWt (octave_idx_type nc, octave_idx_type il, octave_idx_type ir, double l, double rr) + CollocWt (octave_idx_type nc, octave_idx_type il, octave_idx_type ir, + double l, double rr) : n (nc), inc_left (il), inc_right (ir), lb (l), rb (rr), - Alpha (0.0), Beta (0.0), r (), q (), A (), B (), initialized (0) { } + Alpha (0.0), Beta (0.0), r (), q (), A (), B (), initialized (false) { } - CollocWt (octave_idx_type nc, double a, double b, octave_idx_type il, octave_idx_type ir) + CollocWt (octave_idx_type nc, double a, double b, octave_idx_type il, + octave_idx_type ir) : n (nc), inc_left (il), inc_right (ir), lb (0.0), rb (1.0), - Alpha (a), Beta (b), initialized (0) { } + Alpha (a), Beta (b), initialized (false) { } - CollocWt (octave_idx_type nc, double a, double b, octave_idx_type il, octave_idx_type ir, + CollocWt (octave_idx_type nc, double a, double b, octave_idx_type il, + octave_idx_type ir, double ll, double rr) : n (nc), inc_left (il), inc_right (ir), lb (ll), rb (rr), - Alpha (a), Beta (b), r (), q (), A (), B (), initialized (0) { } + Alpha (a), Beta (b), r (), q (), A (), B (), initialized (false) { } CollocWt (const CollocWt& a) : n (a.n), inc_left (a.inc_left), inc_right (a.inc_right), @@ -85,21 +88,21 @@ CollocWt& resize (octave_idx_type nc) { n = nc; - initialized = 0; + initialized = false; return *this; } CollocWt& add_left (void) { inc_left = 1; - initialized = 0; + initialized = false; return *this; } CollocWt& delete_left (void) { inc_left = 0; - initialized = 0; + initialized = false; return *this; } @@ -108,14 +111,14 @@ CollocWt& add_right (void) { inc_right = 1; - initialized = 0; + initialized = false; return *this; } CollocWt& delete_right (void) { inc_right = 0; - initialized = 0; + initialized = false; return *this; } @@ -124,14 +127,14 @@ CollocWt& set_alpha (double val) { Alpha = val; - initialized = 0; + initialized = false; return *this; } CollocWt& set_beta (double val) { Beta = val; - initialized = 0; + initialized = false; return *this; } @@ -178,7 +181,7 @@ Matrix A; Matrix B; - int initialized; + bool initialized; void init (void);