Initial commit from Polytechnique Montreal

1 files changed, 157 insertions, 0 deletions

diff --git a/Utilib/src/ALPRTB.f b/Utilib/src/ALPRTB.f
new file mode 100644
index 0000000..068066c
--- /dev/null
+++ b/Utilib/src/ALPRTB.f
@@ -0,0 +1,157 @@
+*DECK ALPRTB
+      SUBROUTINE ALPRTB(NOR,IINI,DEMT,IER,WEIGHT,BASEPT)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* compute a probability table preserving 2*NOR moments of a function
+* using the modified Ribon approach.
+*
+*Copyright:
+* Copyright (C) 1993 Ecole Polytechnique de Montreal
+* This library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public
+* License as published by the Free Software Foundation; either
+* version 2.1 of the License, or (at your option) any later version
+*
+*Author(s): A. Hebert
+*
+*Parameters: input
+* NOR     the number of moments to preserve is 2*NOR.
+* IINI    minimum order of the moment we want to preserve. we must
+*         have 2-2*NOR <= IINI <= 0 (order 0 and 1 moments are always
+*         preserved).
+* DEMT    moments.
+*
+*Parameters: output
+* IER     error flag (=0/=1 success/failure of the algorithm).
+* WEIGHT  weights of the probability table.
+* BASEPT  base points of the probability table.
+*
+*-----------------------------------------------------------------------
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER NOR,IINI,IER
+      DOUBLE PRECISION DEMT(IINI:2*NOR+IINI-1)
+      REAL WEIGHT(NOR),BASEPT(NOR)
+*----
+*  LOCAL VARIABLES
+*----
+      PARAMETER (MAXNOR=20)
+      DOUBLE PRECISION DS(MAXNOR+1,MAXNOR+1),DDA(0:MAXNOR),DD,DSIGX
+      COMPLEX*16 ROOTS(MAXNOR),CCC,DCC,XCC
+      COMPLEX CGAR
+      LOGICAL LFAIL
+*
+      IF(NOR.GT.MAXNOR) CALL XABORT('ALPRTB: STORAGE OVERFLOW.')
+      IF(NOR.LE.0) CALL XABORT('ALPRTB: NEGATIVE OR ZERO VALUE OF NOR.')
+      IF((2-2*NOR.GT.IINI).OR.(IINI.GT.0)) CALL XABORT('ALPRTB: INCONSI'
+     1 //'STENT VALUE OF IINI.')
+*
+* BUILD THE MATRIX.
+      DO 15 IOR=1,NOR
+      DS(IOR,NOR+1)=-DEMT(NOR+IOR+IINI-1)
+      DO 10 JOR=1,IOR
+      DS(IOR,JOR)=DEMT(IOR+JOR+IINI-2)
+      DS(JOR,IOR)=DEMT(IOR+JOR+IINI-2)
+   10 CONTINUE
+   15 CONTINUE
+*
+* L-D-L(T) FACTORIZATION OF THE MATRIX.
+      DO 40 I=1,NOR
+      DO 30 J=1,I-1
+      DS(J,I)=DS(I,J)
+      DO 20 K=1,J-1
+      DS(J,I)=DS(J,I)-DS(K,I)*DS(J,K)
+   20 CONTINUE
+      DS(I,J)=DS(J,I)*DS(J,J)
+      DS(I,I)=DS(I,I)-DS(J,I)*DS(I,J)
+   30 CONTINUE
+      IF(DS(I,I).EQ.0.D0) THEN
+         IER=1
+         RETURN
+      ENDIF
+      DS(I,I)=1.D0/DS(I,I)
+   40 CONTINUE
+*
+* SOLUTION OF THE FACTORIZED SYSTEM TO OBTAIN THE DENOMINATOR OF THE
+* PADE APPROXIMATION.
+      DO 55 I=1,NOR
+      DO 50 K=1,I-1
+      DS(I,NOR+1)=DS(I,NOR+1)-DS(I,K)*DS(K,NOR+1)
+   50 CONTINUE
+   55 CONTINUE
+      DO 60 I=1,NOR
+      DS(I,NOR+1)=DS(I,NOR+1)*DS(I,I)
+   60 CONTINUE
+      DO 71 I=NOR,1,-1
+      DO 70 K=I+1,NOR
+      DS(I,NOR+1)=DS(I,NOR+1)-DS(K,I)*DS(K,NOR+1)
+   70 CONTINUE
+   71 CONTINUE
+      DS(NOR+1,NOR+1)=1.0D0
+*
+* COMPUTE THE BASE POINTS AS THE ROOTS OF THE DENOMINATOR.
+      CALL ALROOT(DS(1,NOR+1),NOR,ROOTS,LFAIL)
+      IF(LFAIL) CALL XABORT('ALPRTB: POLYNOMIAL ROOT FINDING FAILURE.')
+      DO 80 I=1,NOR
+*
+*     NEWTON IMPROVEMENT OF THE ROOTS.
+      CCC=0.0D0
+      XCC=1.0D0
+      DO 74 J=0,NOR
+      CCC=CCC+DS(J+1,NOR+1)*XCC
+      XCC=XCC*ROOTS(I)
+   74 CONTINUE
+      DCC=0.0D0
+      XCC=1.0D0
+      DO 75 J=1,NOR
+      DCC=DCC+DS(J+1,NOR+1)*XCC*REAL(J)
+      XCC=XCC*ROOTS(I)
+   75 CONTINUE
+      ROOTS(I)=ROOTS(I)-CCC/DCC
+*
+      CGAR=CMPLX(ROOTS(I))
+      IF(ABS(AIMAG(CGAR)).GT.1.0E-4*ABS(REAL(CGAR))) THEN
+         IER=1
+         RETURN
+      ELSE
+         BASEPT(I)=REAL(CMPLX(ROOTS(I)))
+      ENDIF
+   80 CONTINUE
+*
+* COMPUTE THE WEIGHTS.
+      DO 130 I=1,NOR
+      DSIGX=DBLE(ROOTS(I))
+      DDA(0)=1.0D0
+      J0=0
+      DO 100 J=1,NOR
+      IF(J.EQ.I) GO TO 100
+      J0=J0+1
+      DDA(J0)=DDA(J0-1)
+      DO 90 K=1,J0-1
+      DDA(J0-K)=DDA(J0-K-1)-DDA(J0-K)*DBLE(ROOTS(J))
+   90 CONTINUE
+      DDA(0)=-DDA(0)*DBLE(ROOTS(J))
+  100 CONTINUE
+      DD=0.0D0
+      DO 110 J=0,NOR-1
+      DD=DD+DDA(J)*DEMT((IINI-1)/2+J)
+  110 CONTINUE
+      DO 120 J=1,NOR
+      IF(J.NE.I) DD=DD/(DBLE(ROOTS(J))-DSIGX)
+  120 CONTINUE
+      WEIGHT(I)=REAL(((-1.0D0)**(NOR-1))*DD*DSIGX**((1-IINI)/2))
+  130 CONTINUE
+*
+* TEST THE CONSISTENCY OF THE SOLUTION.
+      DO 140 I=1,NOR
+      IF((WEIGHT(I).LE.0.0).OR.(BASEPT(I).LE.0.0)) THEN
+         IER=1
+         RETURN
+      ENDIF
+  140 CONTINUE
+      IER=0
+      RETURN
+      END