diff options
Diffstat (limited to 'Utilib/src/ALPLSF.f')
| -rw-r--r-- | Utilib/src/ALPLSF.f | 282 |
1 files changed, 282 insertions, 0 deletions
diff --git a/Utilib/src/ALPLSF.f b/Utilib/src/ALPLSF.f new file mode 100644 index 0000000..53dc616 --- /dev/null +++ b/Utilib/src/ALPLSF.f @@ -0,0 +1,282 @@ +*DECK ALPLSF + SUBROUTINE ALPLSF(IMETH,N,X,Y,EPSRID,LREAL,NOR,A,B,PREC) +* +*----------------------------------------------------------------------- +* +*Purpose: +* compute the polynomial coefficients of a Pade approximation using a +* direct least square procedure. +* +*Copyright: +* Copyright (C) 1996 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 +* IMETH type of algorithm (=1: use QR; =2: use SVD; =3: use NNLS). +* N number of collocation points. +* X abscissa of the collocation points. +* Y ordinates of the collocation points. +* EPSRID epsilon used in polynomial simplification. +* LREAL selection flag (=.true. to get rid of complex roots). +* +*Parameters: output +* NOR order of the polynomials. +* A polynomial coefficients of the numerator of the Pade +* approximation. A(0) is the constant term. +* B polynomial coefficients of the denominator of the Pade +* approximation. B(0) is the constant term. +* DOUBLE PRECISION A(0:NOR),B(0:NOR) +* PREC accuracy of the fit. +* +*----------------------------------------------------------------------- +* +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER IMETH,N,NOR + REAL X(N),Y(N),PREC + DOUBLE PRECISION EPSRID,A(0:(N-1)/2),B(0:(N-1)/2) + LOGICAL LREAL +*---- +* LOCAL VARIABLES +*---- + PARAMETER (MAXNOR=10,MAXPTS=99) + DOUBLE PRECISION BB(MAXPTS),GAR,PARAM(MAXPTS),AGAR(0:MAXNOR), + 1 BGAR(0:MAXNOR),CGAR(0:MAXNOR+1),W(MAXPTS),RV1(MAXPTS),SGN, + 2 GAROLD,YAPPR,RNORM,RMAX + COMPLEX*16 SIGX0(MAXNOR+1),SIGXW(MAXNOR+1),DDAGAR(0:MAXNOR), + 1 DDBGAR(0:MAXNOR),WEIGH(MAXNOR+1),CC,DD,CCC,XCC,DCC + LOGICAL LFAIL + DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: AA,V +*---- +* SCRATCH STORAGE ALLOCATION +*---- + ALLOCATE(AA(MAXPTS,MAXPTS),V(MAXPTS,MAXPTS)) +* + IF(N.GT.MAXPTS) CALL XABORT('ALPLSF: INSUFFICIENT MAXPTS.') +* +* NOR=0 CASE. + GAROLD=0.0 + NOR=0 + PREC=0.0 + A(0)=DBLE(Y(N)) + B(0)=1.0D0 + DO 10 I=1,N-1 + PREC=MAX(PREC,ABS(Y(N)/Y(I)-1.0)) + 10 CONTINUE +* + RMAX=1.0D10 + DO 210 IINOR=1,MIN((N-1)/2,MAXNOR) + INOR=IINOR + IF(X(N).GE.0.99E10) THEN + DO 41 I=1,N-1 + GAR=1.0D0 + IOF=0 + DO 20 J=0,INOR-1 + IOF=IOF+1 + AA(I,IOF)=GAR + GAR=GAR*X(I) + 20 CONTINUE + GAR=1.0D0 + DO 30 J=0,INOR-1 + IOF=IOF+1 + GAROLD=GAR + AA(I,IOF)=-GAR*Y(I) + GAR=GAR*X(I) + 30 CONTINUE + BB(I)=(Y(I)-Y(N))*X(I) + DO 40 J=1,2*INOR + AA(I,J)=AA(I,J)/GAROLD + 40 CONTINUE + 41 CONTINUE + IF(IMETH.EQ.1) THEN + CALL ALST2F(MAXPTS,N-1,2*INOR,AA,RV1) + CALL ALST2S(MAXPTS,N-1,2*INOR,AA,RV1,BB,PARAM) + ELSE IF(IMETH.EQ.2) THEN + CALL ALSVDF(AA,N-1,2*INOR,MAXPTS,MAXPTS,W,V,RV1) + DO 45 J=1,2*INOR + IF(W(J).EQ.0.0D0) CALL XABORT('ALPLSF: SVD FAILURE(1).') + 45 CONTINUE + CALL ALSVDS(AA,W,V,N-1,2*INOR,MAXPTS,MAXPTS,BB,PARAM,RV1) + ELSE IF(IMETH.EQ.3) THEN + CALL ALNNLS(AA,N-1,2*INOR,MAXPTS,MAXPTS,BB,PARAM,RNORM,MODE) + IF(MODE.NE.1) CALL XABORT('ALPLSF: NNLS FAILURE(1).') + IF((INOR.GT.1).AND.(RNORM.GE.0.95D0*RMAX)) GO TO 210 + RMAX=RNORM + ENDIF + DO 50 I=0,INOR-1 + AGAR(I)=PARAM(I+1) + BGAR(I)=PARAM(INOR+1+I) + 50 CONTINUE + AGAR(INOR)=Y(N) + BGAR(INOR)=1.0D0 + ELSE + DO 81 I=1,N + GAR=1.0D0 + IOF=0 + DO 60 J=0,INOR + IOF=IOF+1 + AA(I,IOF)=GAR + GAR=GAR*X(I) + 60 CONTINUE + GAR=1.0D0 + DO 70 J=0,INOR-1 + IOF=IOF+1 + GAROLD=GAR + AA(I,IOF)=-GAR*Y(I) + GAR=GAR*X(I) + 70 CONTINUE + BB(I)=Y(I)*X(I) + DO 80 J=1,2*INOR+1 + AA(I,J)=AA(I,J)/GAROLD + 80 CONTINUE + 81 CONTINUE + IF(IMETH.EQ.1) THEN + CALL ALST2F(MAXPTS,N,2*INOR+1,AA,RV1) + CALL ALST2S(MAXPTS,N,2*INOR+1,AA,RV1,BB,PARAM) + ELSE IF(IMETH.EQ.2) THEN + CALL ALSVDF(AA,N,2*INOR+1,MAXPTS,MAXPTS,W,V,RV1) + DO 85 J=1,2*INOR + IF(W(J).EQ.0.0D0) CALL XABORT('ALPLSF: SVD FAILURE(2).') + 85 CONTINUE + CALL ALSVDS(AA,W,V,N,2*INOR+1,MAXPTS,MAXPTS,BB,PARAM,RV1) + ELSE IF(IMETH.EQ.3) THEN + CALL ALNNLS(AA,N,2*INOR+1,MAXPTS,MAXPTS,BB,PARAM,RNORM,MODE) + IF(MODE.NE.1) CALL XABORT('ALPLSF: NNLS FAILURE(2).') + IF((INOR.GT.1).AND.(RNORM.GE.0.95D0*RMAX)) GO TO 210 + RMAX=RNORM + ENDIF + DO 90 I=0,INOR + AGAR(I)=PARAM(I+1) + IF(I.EQ.INOR) THEN + BGAR(I)=1.0D0 + ELSE + BGAR(I)=PARAM(INOR+2+I) + ENDIF + 90 CONTINUE + ENDIF +* +* POLYNOMIAL SIMPLIFICATION. + DDAGAR(0)=AGAR(INOR) + DDBGAR(0)=BGAR(INOR) + CALL ALROOT(AGAR(0:INOR),INOR,SIGX0,LFAIL) + IF(LFAIL) GO TO 210 + CALL ALROOT(BGAR(0:INOR),INOR,SIGXW,LFAIL) + IF(LFAIL) GO TO 210 + IJINOR=1 + 95 XXX=REAL(ABS(DBLE(SIGXW(IJINOR))-DBLE(SIGX0(IJINOR)))) + IF(XXX.LT.EPSRID*ABS(DBLE(SIGXW(IJINOR)))) THEN + INOR=INOR-1 + DO 100 I=IJINOR,INOR + SIGX0(I)=SIGX0(I+1) + SIGXW(I)=SIGXW(I+1) + 100 CONTINUE + ELSE IF((DBLE(SIGXW(IJINOR)).GT.EPSRID).AND. + > (DIMAG(SIGXW(IJINOR)).EQ.0.0).AND. + > (IMETH.EQ.3)) THEN + CALL XABORT('ALPLSF: NNLS FAILURE(3).') + ELSE IF((DBLE(SIGXW(IJINOR)).GT.0.1*EPSRID).AND. + > (DIMAG(SIGXW(IJINOR)).EQ.0.0)) THEN + GO TO 210 + ELSE + IJINOR=IJINOR+1 + ENDIF + IF(IJINOR.LE.INOR) GO TO 95 + IF(INOR.LT.0) CALL XABORT('ALPLSF: ALGORITHM FAILURE.') + DO 120 I=1,INOR + DDAGAR(I)=DDAGAR(I-1) + DDBGAR(I)=DDBGAR(I-1) + DO 110 J=I-1,1,-1 + DDAGAR(J)=DDAGAR(J-1)-DDAGAR(J)*SIGX0(I) + DDBGAR(J)=DDBGAR(J-1)-DDBGAR(J)*SIGXW(I) + 110 CONTINUE + DDAGAR(0)=-DDAGAR(0)*SIGX0(I) + DDBGAR(0)=-DDBGAR(0)*SIGXW(I) + 120 CONTINUE + DO 130 I=0,INOR + AGAR(I)=DBLE(DDAGAR(I))/DBLE(DDBGAR(INOR)) + BGAR(I)=DBLE(DDBGAR(I))/DBLE(DDBGAR(INOR)) + IF(AGAR(I).LE.0.0D0) GO TO 210 + IF(BGAR(I).LE.0.0D0) GO TO 210 + 130 CONTINUE + SGN=1.0D0 + CGAR(0)=AGAR(0) + DO 135 I=2,INOR+1 + SGN=-SGN + CGAR(I-1)=SGN*(BGAR(I-2)+AGAR(I-1)) + 135 CONTINUE + CGAR(INOR+1)=-SGN + CALL ALROOT(CGAR(0:INOR+1),INOR+1,SIGX0,LFAIL) + IF(LFAIL) GO TO 210 +* +* NEWTON IMPROVEMENT OF THE ROOTS. + DO 138 I=1,INOR+1 + CCC=0.0D0 + XCC=1.0D0 + DO 136 J=0,INOR+1 + CCC=CCC+CGAR(J)*XCC + XCC=XCC*SIGX0(I) + 136 CONTINUE + DCC=0.0D0 + XCC=1.0D0 + DO 137 J=1,INOR+1 + DCC=DCC+CGAR(J)*XCC*REAL(J) + XCC=XCC*SIGX0(I) + 137 CONTINUE + SIGX0(I)=SIGX0(I)-CCC/DCC + 138 CONTINUE +* + IF(LREAL) THEN + DO 140 I=1,INOR+1 + IF(DBLE(SIGX0(I)).LT.1.0E-10) GO TO 210 + IF(DIMAG(SIGX0(I)).NE.0.0) GO TO 210 + 140 CONTINUE + ENDIF +* +* COMPUTE THE WEIGHTS. + DO 170 I=1,INOR+1 + CC=(1.0D0,0.0D0) + DD=0.0D0 + DO 150 JNOR=0,INOR + DD=DD+BGAR(JNOR)*CC + CC=-CC*SIGX0(I) + 150 CONTINUE + DO 160 J=1,INOR+1 + IF(J.NE.I) DD=DD/(SIGX0(J)-SIGX0(I)) + 160 CONTINUE + WEIGH(I)=DD + 170 CONTINUE +* +* TEST THE ACCURACY OF THE PADE APPROXIMATION. + PREC1=0.0 + DO 190 I=1,N + CC=0.0D0 + DD=0.0D0 + DO 180 JNOR=1,INOR+1 + CC=CC+WEIGH(JNOR)/(SIGX0(JNOR)+X(I)) + DD=DD+WEIGH(JNOR)*SIGX0(JNOR)/(SIGX0(JNOR)+X(I)) + 180 CONTINUE + YAPPR=DBLE(DD/CC) + PREC1=MAX(PREC1,ABS(REAL(YAPPR)/Y(I)-1.0)) + 190 CONTINUE +* + IF(PREC1.LT.0.95*PREC) THEN + NOR=INOR + PREC=PREC1 + DO 200 I=0,NOR + A(I)=AGAR(I) + B(I)=BGAR(I) + 200 CONTINUE + ENDIF + 210 CONTINUE +*---- +* SCRATCH STORAGE DEALLOCATION +*---- + DEALLOCATE(V,AA) + RETURN + END |