Differences between revisions 185 and 186
Revision 185 as of 2011-06-25 22:32:31
Size: 6832
Editor: 76
Comment:
Revision 186 as of 2011-08-07 15:39:55
Size: 13184
Editor: 221
Comment:
Deletions are marked like this. Additions are marked like this.
Line 1: Line 1:
= Welcome to GCC Wiki =
This page contains information about the [[http://gcc.gnu.org|GNU Compiler Collection]]. Please read HowToUseWiki if you need help editing this WikiWikiWeb.
C 19-TH February, 1992
C
C INPUT R* AND V INSTEAD MU AND BETA
C TRY TO FIND THE CRITICAL GEOMETRIC RATIO FOR RUBBER
C
C EIGENVALUES SURFACE OF ANTISYMMETRIC DIFFUSE MODES
C ************ FOR THE EI SUBREGIME ***************
C
C FOR INCOMPRESSIBLE ISOTROPIC PRESSURE-INSENSITIVE MATERIAL (RUBBER)
C UNDER NO CONFINING STRESSES
C
C TRY TO COMPARE WITH BEATTY'S (1976) EXPERIMENTAL RESULTS
C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
C FOR PLOTTING ALL THOSE ELLIPTIC EIGENVALUES MODE
C WITH CHANGING SZZ/(2.*GL) RATIO AND GAMMA (WAVELENGTH)
C
C DEVELOPED BY
C K. T. CHAU
C
C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
 PROGRAM QAMPAPER
 IMPLICIT REAL*8 (A-H,O-Z)
 DIMENSION XB1(10),XB2(10)
 EXTERNAL FEIS,FEIA,ZBRENT
  CALL CTIME(0)
C
C HAVE TO INPUT FIVE MATERIAL PARAMETERS
C AND ONE PARAMETER IS CHANGING WITH CURVES
C
C READ(5,10001)RG,RGK,PR,RR,SRE
C READ(5,10002)REG
C10001 FORMAT(5F10.0)
C10002 FORMAT(F10.0)
C
C WE HAVE
 RG =1.0
 RGK=0.0
 PR =0.5
 RR =1.0
 SRE=0.0
 REG=3.0
C IN THIS PROGRAM
C
C DEFINE THE FOLLOWING RATIO :
C
C PR = POISSON RATIO
C R = PRESSURE-SENSITIVE COEFFICIENT
C RG = GT/GL
C RGK = GL/K
C SRE = SRR/E
C SZE = SZZ/E (THIS THE UNKNOWN TO BE SOLVED)
C RVV = R/2V (NON-NORMALITY PARAMETER)
C WHERE
C GL = LONGITUDINAL SHEAR MODULUS
C GT = TRANSVERSE SHEAR MODULUS
C E = INSTANTANEOUS TANGENT MODULUS
C K = IN-PLANE BULK MODULUS
C RR = PRESSURE-SENSITIVITY FACTOR
C PR = EF
C EFFCTIVE POISSON'S RATIO (V)
C
 DO 100 KK1=1,2
C
C KK1=1 FOR SYMMETRIC MODE
C KK2=2 FOR ANTISYMMETRIC MODE
C
 WRITE(6,20008)RG,RGK,PR,RR,REG,SRE
C
20008 FORMAT(5X,' APPLIED COEFFICENT OF THE PROBLEM'/
     * 5X,' GT/GL (IN-PLANE BULK MODULUS) =',E12.3/
     * 5X,' GL/K (INSTANTANEOUS E) =',E12.3/
     * 5X,' PR (POISSON RATIO ) =',E12.3/
     * 5X,' RR (PRESSURE-SENSITIVITY) =',E12.3/
     * 5X,' E/GL (BETA INRR MODEL) =',E12.3/
     * 5X,' SRR/E (RADIAL STRESS ) =',E12.3//)
C
C
  TOL=1.E-6
C
C GAMMA = N(PI)A/L
C WHERE
C N = WAVE NUMBER
C A = RADIUS OF THE CYLINDER
C L = LENGTH OF THE CYLINDER
C (GAMMA WILL BE THE X-AXIS OF THE PLOT)
C
 GAMMA=50.0
  DG=GAMMA/100.
C
C
    WRITE(6,2011)REG
2011 FORMAT(1X,' E/GL =',E20.6//)
C
 DO 200 KK2=1,1
C
C KK2 =1 FOR S<0
C KK2 =2 FOR S>0
C
C
 WRITE(6,292)
 WRITE(0,292)
292 FORMAT(1X,' GAM SZZ/2GL',//)
   DO 221 IIB=1,101
   GAM=(IIB-1)*DG
C
C ELLIPTIC COMPLEX Regime
C
 X1=0.00
   X2=-0.99
  IF(KK2.EQ.2) X2=0.99
   N=10
C
C NB = NO. OF ROOTS WANTED
C N= NO. DIVISION INPUT
C
   NB=2
 IF(KK1.EQ.1) THEN
  CALL ZBRAK(FEIS,X1,X2,GAM,XB1,XB2,N,NB)
 ELSE
  CALL ZBRAK(FEIA,X1,X2,GAM,XB1,XB2,N,NB)
 ENDIF
C
C
  IF(NB.EQ.0) THEN
 WRITE(0,12011)GAM
12011 FORMAT(1X,'THERE IS NO ROOT (EC REGIME) AT GAM=',F20.5)
  GO TO 221
 ENDIF
C XACC=1.E-6
 DO 776 LL1=1,NB
  XX1=XB1(LL1)
  XX2=XB2(LL1)
 IF(KK1.EQ.1) THEN
  SZZ2=ZBRENT(FEIS,XX1,XX2,GAM)
 ELSE
  SZZ2=ZBRENT(FEIA,XX1,XX2,GAM)
 ENDIF
C
 WRITE(6,21111)GAM,SZZ2
 WRITE(0,21111)GAM,SZZ2
776 CONTINUE
21111 FORMAT(2E20.6)
C
221 CONTINUE
C
C
C Elliptic-Imaginary
C
C
C--------------------------------------------
C END OF THE EGL (VARYING) LOOP
C---------------------------------------
200 CONTINUE
100 CONTINUE
999 CONTINUE
 CALL CTIME (-1)
 STOP
 END
C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C DEVELOPED BY
C K. T. CHAU
C REF: NUMERICAL RECIPES,PP.252
C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C Using Brent's method, find the root of a function FUNC knowm to
c lie between X1 and X2. The root, returned as ZBREET, will be
c refined until its accuracy is TOL.
C
 FUNCTION ZBRENT(FUNCC,X3,X4,GAM)
 IMPLICIT REAL*8(A-H,O-Z)
 EXTERNAL FUNCC
 ITMAX=250
 TOL=1.E-6
 EPS=1.E-16
 A=X3
 B=X4
C
 FA=FUNCC(A,GAM)
 FB=FUNCC(B,GAM)
CC
C
 IF (FB*FA.GT.0.) PAUSE 'ROOT MUST BE BRACKETED FOR ZBRENT!'
 FC=FB
 DO 11 ITER=1,ITMAX
         IF (FB*FC.GT.0.0) THEN
C
C RENAME A, B AND C AND ADJUST BOUNDING INTERVAL D.
C
 C=A
 FC=FA
 D=B-A
 E=D
 ENDIF
 IF (DABS(FC).LT.DABS(FB)) THEN
 A=B
 B=C
 C=A
 FA=FB
 FB=FC
 FC=FA
 ENDIF
CC
C CONVERGENCE CHECK
C
 TOL1=2.*EPS*DABS(B)+0.5*TOL
 XM=0.5*(C-B)
 IF (DABS(XM).LE.TOL1.OR.FB.EQ.0.0) THEN
 ZBRENT=B
 RETURN
 ENDIF
 IF (DABS(E).GE.TOL1.AND.DABS(FA).GT.DABS(FB)) THEN
CC ATTEMPT INVERSE QUADRATIC INTERPOLATION
C
 S=FB/FA
 IF (A.EQ.C) THEN
 P=2.*XM*S
 Q=1.-S
 ELSE
 Q=FA/FC
 R=FB/FC
 P=S*(2.*XM*Q*(Q-R)-(B-A)*(R-1.))
 Q=(Q-1.)*(R-1.)*(S-1.)
 ENDIF
C
C CHECK WHETHER IN BOUNDS
C
 IF (P.GT.0.0) Q=-Q
 P=DABS(P)
 IF (2.*P.LT.DMIN1(3.*XM*Q-DABS(TOL1*Q),DABS(E*Q))) THEN
C
C ACCEPT INTERPOLATION
C
 E=D
 D=P/Q
 ELSE
C
C INTERPOLATION FAILED, USE BISECTION.
C
 D=XM
 E=D
 ENDIF
 ELSE
C
C BOUNDS DECREASES TOO SLOWLY, USE BISECTION.
C
 D=XM
 E=D
 ENDIF
C
C MOVE LAST BEST GUESS TO A
C
 A=B
 FA=FB
 IF(DABS(D).GT.TOL1) THEN
C
C EVALUATE NEW TRIAL ROOT
C
 B=B+D
 ELSE
 B=B+DSIGN(TOL1,XM)
 ENDIF
C
 FB=FUNCC(B,GAM)
C
11 CONTINUE
 PAUSE 'ZBRENT EXCEEDING MAXIMUM ITERATIONS'
 ZBRENT=B
   RETURN
 END
C
C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C SUBROUTINE TO BRACKET THE INTERVAL FROM X1 TO X2
C
 SUBROUTINE ZBRAK(FX,X1,X2,GAM,XBB1,XBB2,N,NB)
C
C GIVEN A FUNCTION FX DEFINED ON THE INTERVAL FROM X1 TO X2
C SUBDIVIDE THE INTERVAL INTO N EQUALLY SPACED SEGMENTS, AND
C SEARCH FOR ZERO CROSSINGS OF THE FUNCTION. NB IS INPUT AS
C THE MAXIMUM NUMBER OF ROOTS SOUGHT, AND IS RESET TO THE
C NUMBER OF BRACKETING PAIRS XB1, XB2 THAT ARE FOUND.
C
C
 IMPLICIT REAL*8 (A-H,O-Z)
 DIMENSION XBB1(10),XBB2(10)
 EXTERNAL FX
 NBB=NB
 NB=0
 X=X1
 DX=(X2-X1)/N
C
 FP=FX(X,GAM)
C
C LOOP OVER ALL INTERVALS
C
 DO 11 I=1,N
         X=X+DX
C
  FC=FX(X,GAM)
C
 IF(FC*FP.LT.0.0) THEN
C IF SIGN CHANGE OCCURS THEN RECORD VALUES FOR THE BOUNDS
C
 NB=NB+1
 XBB1(NB)=X-DX
  XBB2(NB)=X
 ENDIF
 FP=FC
 IF (NBB.EQ.NB) RETURN
11 CONTINUE
 RETURN
 END
C---|----1----|----2----|----3----|----4----|----5----|----6----|----712
C
 FUNCTION FEIA(SZZ,GAM)
 IMPLICIT REAL*8 (A-H,O-Z)
  EXTERNAL BESSI0,BESSI1,BESSI
CC
C NEW SUBROUTINE FOR ANTI-SYMMETRIC MODES
C
 SS=SZZ
 VA=(1.-SS)
 VB=2.
 VC=(1.+SS)
C V3 FOR THE SECOND DIFFERENTIAL EQUATION
 U3=DSQRT(1.+SS)
 U32=U3*U3
C
C Elliptic-Complex Regime
C
C IF((VB*VB-4.*VA*VC).LT.0.0.AND.(VA*VC).GT.0.0) THEN
C IF(VA.EQ.0.0) PAUSE 'A=0 IN FECA'
C
C CHECK POSITIVENESS OF P AND Q
C
C CALCULATION OF P AND Q
 PP=DSQRT((VB+DSQRT(VB*VB-4.*VA*VC))/(2.*VA))
 QQ=DSQRT((VB-DSQRT(VB*VB-4.*VA*VC))/(2.*VA))
 PP2=PP*PP
 QQ2=QQ*QQ
 GP=GAM*PP
 GQ=GAM*QQ
C
 BIP0=BESSI0(GP)
 BIP1=BESSI1(GP)
 BIP2=BESSI(2,GP)
 BIP3=BESSI(3,GP)
 BIQ0=BESSI0(GQ)
 BIQ1=BESSI1(GQ)
 BIQ2=BESSI(2,GQ)
 BIQ3=BESSI(3,GQ)
C
 ARG=GAM*U3
 BI0 =BESSI0(ARG)
 BI1 =BESSI1(ARG)
 BI2 =BESSI(2,ARG)
 BI3 =BESSI(3,ARG)
C---|----1----|----2----|----3----|----4----|----5----|----6----|----712
 Y1 = -(1.+SS)
 Y2 = Y1
 Y3 = 2.0
 Y4 = -(1.-SS)
 Y5 = (1.-SS)
 Y6 = Y5
 Y7 = 1.0
 Y8 = Y7
C
 IF(GAM.EQ.0.0) GO TO 777
 WW1 = -4.*U32*Y3*BI2/ARG
 WW2 = U3*Y6*BI1/ARG
 WW3 = U32/4.*(2.*Y7/ARG*(BI0+BI2-2.*BI1/ARG)-Y8*(3.*BI1+BI3))
C
 VPP=((4.*Y1+3.*Y3)*PP2-4.*Y2)*BIP1+PP2*Y3*BIP3
 VQP=0.5*PP*(Y4*PP2-Y5)*(BIP0+BIP2)
 VRP=PP2*Y3*BIP2/GP
 VPQ=((4.*Y1+3.*Y3)*QQ2-4.*Y2)*BIQ1+QQ2*Y3*BIQ3
 VQQ=0.5*QQ*(Y4*QQ2-Y5)*(BIQ0+BIQ2)
 VRQ=QQ2*Y3*BIQ2/GQ
C
C---|----1----|----2----|----3----|----4----|----5----|----6----|----712
 GO TO 778
777 FEIA=U32*Y3*Y3*Y5+4.*Y2*Y3*Y6+U32*(Y7-3.*Y8)*(Y5*(4.*Y1+3.*Y3)
     * -4.*Y2*Y4)
 GO TO 999
C ELSE
778 FEIA=WW1*(VQP*VRQ-VQQ*VRP)+WW2*(VPQ*VRP-VPP*VRQ)
     * +WW3*(VPP*VQQ-VPQ*VQP)
999 RETURN
 END
C-------------------------------------------------------------------------
C EIGENVALUE EQUATION FOR SYMMETRIC GEOMETRIC DIFFUSE MODE
C-------------------------------------------------------------------------
 FUNCTION FEIS(SZZ,GAM)
 IMPLICIT REAL*8 (A-H,O-Z)
 EXTERNAL BESSI1,BESSI0,BESSI
C
 SS=SZZ
 VA=(1.-SS)
 VB=2.
 VC=(1.+SS)
C
C Elliptic-Imaginary Subregime
C
 PP=DSQRT((VB+DSQRT(VB*VB-4.*VA*VC))/(2.*VA))
 QQ=DSQRT((VB-DSQRT(VB*VB-4.*VA*VC))/(2.*VA))
c PP=DSQRT((2.*DSQRT(VA*VC)-VB)/(4.*VA))
c QQ=DSQRT((2.*DSQRT(VA*VC)+VB)/(4.*VA))
 PP2=PP*PP
 QQ2=QQ*QQ
 GP=GAM*PP
 GQ=GAM*QQ
C
 BIP0=BESSI0(GP)
 BIP1=BESSI1(GP)
 BIP2=BESSI(2,GP)
 BIQ0=BESSI0(GQ)
 BIQ1=BESSI1(GQ)
 BIQ2=BESSI(2,GQ)
C---|----1----|----2----|----3----|----4----|----5----|----6----|----712
 Y1 = -(1.+SS)
 Y2 = Y1
 Y3 = 2.0
 Y4 = -(1.-SS)
 Y5 = (1.-SS)
 Y6 = Y5
 Y7 = 1.0
 Y8 = Y7
C
 VPP=2.*((2.*Y1+Y3)*PP2-2.*Y2)*BIP0+PP2*Y3*2.*BIP2
 VPQ=2.*((2.*Y1+Y3)*QQ2-2.*Y2)*BIQ0+QQ2*Y3*2.*BIQ2
 VQP=2.*PP*(Y4*PP2-Y5)*BIP1
 VQQ=2.*QQ*(Y4*QQ2-Y5)*BIQ1
C
C
 IF(GAM.EQ.0.0) THEN
 FEIS=VC*Y4*(2.*Y1+Y3)+2.*Y2*(VA*Y5-VB*Y4)
 GO TO 999
 ELSE
 FEIS=VPP*VQQ-VPQ*VQP
 ENDIF
999 RETURN
 END
C
C SUBROUTINE FOR TIME CALCULATION FOR THE RUNNING TIME
C ****************************************************
C
  INTERFACE TO SUBROUTINE TIME (N,STR)
 CHARACTER*10 STR [NEAR,REFERENCE]
 INTEGER*2 N [VALUE]
 END
$LARGE
 SUBROUTINE CTIME(N)
        IMPLICIT REAL*8 (A-H,O-Z)
 CHARACTER*1 C1,C2,C4,C5,C7,C8
        CHARACTER*10 CLOCK
        COMMON/CPUT/CT,PT
 CALL TIME (10,CLOCK)
 C1=CLOCK(1:1)
 C2=CLOCK(2:2)
 C4=CLOCK(4:4)
 C5=CLOCK(5:5)
 C7=CLOCK(7:7)
 C8=CLOCK(8:8)
 HOUR=(ICHAR(C1)-48)*10+ICHAR(C2)-48
 MINU=(ICHAR(C4)-48)*10+ICHAR(C5)-48
 SECO=(ICHAR(C7)-48)*10+ICHAR(C8)-48
 A=HOUR*3600.+MINU*60.+SECO
        IF ( N.EQ.0) GO TO 100
        IF ( N.LT.0) GO TO 200
        IF ( N.GT.0) GO TO 300
300 CT = A
        IF (CT.LT.PT) CT = CT + 3600.*24.
        T = CT - PT
        WRITE (6,1000) T
        WRITE (0,1000) T
1000 FORMAT (41X, ' CURRENT TIME (SEC.)>>>>' ,F10.0,/)
        RETURN
100 PT = A
        ST = 0.0
        WRITE (6,2000) ST
        WRITE (0,2000) ST
2000 FORMAT (41X, ' STARTING TIME (SEC.)>>>>',F10.0,/)
 RETURN
200 CT = A
        IF (CT.LT.PT) CT = CT + 3600.*24.
        T = CT - PT
        WRITE (6,3000) T
        WRITE (0,3000) T
3000 FORMAT (41X, ' ENDING TIME (SEC.)>>>>',F10.0,/)
        RETURN
        END
C---|----1----|----2----|----3----|----4----|----5----|----6----|----712
 FUNCTION BESSI0(X)
 IMPLICIT REAL*8 (A-H,O-Z)
 DATA P1,P2,P3,P4,P5,P6,P7/1.D0,3.5156229D0,3.0899424D0,
     * 1.2067492D0,0.2659732D0,0.360768D-1,0.45813D-2/
        DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9/0.39894228D0,0.1328592D-1,
     * 0.225319D-2,-0.157565D-2,0.916281D-2,-0.2057706D-1,
     * 0.2635537D-1,-0.1647633D-1,0.392377D-2/
        IF(DABS(X).LT.3.75) THEN
                Y=(X/3.75)**2.
                BESSI0=P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7)))))
        ELSE
                AX=DABS(X)
                Y=3.75/AX
                BESSI0=(DEXP(AX)/DSQRT(AX))*(Q1+Y*(Q2+Y*(Q3+Y*(Q4
     * +Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9))))))))
        ENDIF
        RETURN
        END
C----------
 FUNCTION BESSI1(X)
 IMPLICIT REAL*8 (A-H,O-Z)
C---|----1----|----2----|----3----|----4----|----5----|----6----|----712
 DATA P1,P2,P3,P4,P5,P6,P7/0.5D0,0.87890594D0,0.51498869D0,
     * 0.15084934D0,0.2658733D-1,0.301532D-2,0.32411D-3/
 DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9/0.39894228D0,-0.3988024D-1,
     * -0.362018D-2,0.163801D-2,-0.1031555D-1,0.2282967D-1,
     * -0.2895312D-1,0.1787654D-1,-0.420059D-2/
        IF(DABS(X).LT.3.75) THEN
                Y=(X/3.75)**2.
Line 4: Line 523:
~+'''This wiki is not for random discussion of GCC, nor for asking questions. It is here to provide information. If you have questions, please use the mailing list. Do not add questions to these pages! gcc-help@gcc.gnu.org is a good mailing list for questions about GCC. '''+~
||<tablewidth="100%">There is an [[GCConIRC|IRC channel dedicated to GCC development]] at {{{irc.oftc.net/#gcc}}} <<BR>> '''NOTE!''' This channel is exclusively dedicated to the '''development''' of GCC. Questions regarding usage of GCC are not appropriate. ||
||There is an IRC channel dedicated to helping users with GCC at {{{irc.freenode.org/#gcc}}}. Please use this channel for questions regarding GCC usage and general GCC discussion. ||




== GCC Info ==
 * [[DevelopmentSchedule|Development schedule]]
 * [[FAQ|Frequently Asked Questions]]
 * [[plugins|GCC plugins]]
 * [[GFortran|GFortran, the Fortran front-end]]
 * [[SummerOfCode|Google's Summer of Code]]
 * [[Libstdc++|libstdc++, the C++ runtime library]]

== Events ==
 * [[GCCGathering2011|GCC Gathering 2011, London, UK for the weekend of 17-Jun-2011]]
 * [[Release Manager Q&A]]
 * [[http://www.gccsummit.org/2010/|GCC Summit 2010, Oct. 25-27, Ottawa, Canada]]
 * [[http://grow2011.inria.fr/|GROW 2011, Apr. 2-3, Chamonix, France, co-located with CGO]]

== Getting Started with GCC Development ==
 * [[GettingStarted|Tutorials, HOWTOs, internal documentation]]
 * [[GCC_Patch_Tracking|GCC Patch Tracking]]
 * [[SvnHelp|SVN Guide for GCC developers]]
 * [[GitMirror|Accessing the GCC sources using Git]]
 * [[Top-Level_Bootstrap|Top-Level Bootstrap]]
 * [[HowToPrepareATestcase|How to prepare a testcase]]
 * [[A_guide_to_testcase_reduction|A guide to testcase reduction]]
 * CompileFarm

== Improving GCC ==
 * BeginnerProjects
 * [[PerformanceTesting|Scripts for testing compile time and memory consumption]]
 * [[Speedup_areas|Speedup areas]]
 * [[general_backend_cleanup|general backend cleanup]]
 * [[Top-Level_Libgcc_Migration|Top-Level Libgcc Migration]]
 * [[Memory_management|Proper GCC Memory Management]]
 * [[Partial_Transitions|Partial Transitions]]
 * BugzillaStats
 * [[openmp|OpenMP support]]
 * [[PDO|Using profile to drive optimizations]]
 * [[Document GCC Internals]]
 * [[gengtype]]

== Current Projects (alphabetical) ==
 * [[Atomic|Atomics]]
 * [[AutoParInGCC|Automatic parallelization]]
 * [[Better_Diagnostics|Better Diagnostics]]
 * [[Better_Uninitialized_Warnings|Better Uninitialized Warnings]]
 * [[Atomic/GCCMM|C++0x Memory Model]]
 * [[Constexpr]]
 * [[http://gcc.gnu.org/projects/cli.html|CLI Back-End and Front-End]]
 * DataflowPorting
 * [[MemRef|Flattening Memory Reference Trees in the GIMPLE IL]]
 * [[functionAdaptation|Function Adaptation]]
 * [[FunctionSpecificOpt|Function Specific Optimization]]
 * [[gcc-in-cxx|gcc-in-cxx - porting gcc to compile as C++]]
 * [[gimplebackend|Gimple Back End]]
 * [[GimpleFrontEnd|Gimple Front End]]
 * [[Graphite]]
 * [[IncrementalCompiler|Incremental compiler]]
 * [[InteractiveCompilationInterface|Interactive Compilation Interface]]
 * [[LightweightIpo|Lightweight IPO -- LIPO]]
 * [[LoopOptTasks|Loop Optimization Related Tasks]]
 * [[NoUndefinedOverflow|Make C undefined overflow semantics explicit in the IL]]
 * [[MiddleEndArrays|Middle End Array Expressions]]
 * [[MiddleEndLispTranslator|Middle End Lisp Translator]] or [[MELT]]
 * [[MilepostGCC|MILEPOST GCC - enabling research on machine-learning based self-tuning compilers]]
 * [[FunctionBehavior|Modeling Function Behavior for more aggressive optimizations across call sites]]
 * [[ModularGCC|Make GCC more modular]]
 * [[SwingModuloScheduling|Modulo Scheduling Related Tasks]]
 * [[OOP|Object Oriented Programming in Fortran]]
 * [[Pass Activity Log]]
 * [[pph|Pre-Parsed Headers]]
 * [[SplitStacks|Split Stacks]]
 * [[SSA Pressure Reduction]]
 * [[PythonFrontEnd|Python Front End]]
 * [[Stdlib Performance Advisor]]
 * [[ThreadSafetyAnnotation|Thread Safety Annotations and Analysis]]
 * [[TransactionalMemory|Transactional Memory]]
 * [[VectorizationTasks|Vectorization Related Tasks]]

== Finished Projects (alphabetical) ==
 * [[Alias_Improvements|Alias Improvements]]
 * [[plugins-branch|Compiler Plugins]]
 * [[tuples|GIMPLE tuples]]
 * [[LinkTimeOptimization|Link Time Optimization]]
 * [[MemoryModel|Tightening GCC's Memory Model]]
 * [[Var_Tracking_Assignments|Var Tracking Assignments for correct debug information]]

== User Information ==
 * [[Visibility|Proper C++ visibility support]]
 * [[Math_Optimization_Flags|Math Optimization Flags]]
 * [[Mudflap_Pointer_Debugging|Mudflap Pointer Debugging]]
 * [[Building_Cross_Toolchains_with_gcc|Building Cross Toolchains with gcc]]
 * [[Software_floating_point|Software floating point]]
 * [[https://twiki.cern.ch/twiki/bin/view/LCG/VILto|Exploiting Link Time Optimization]]

<<Anchor(summitprocs)>>

== GCC Summit Proceedings ==
Proceedings of the [[http://www.gccsummit.org/|GCC Summit]]:

 * [[attachment:2003-GCC-Summit-Proceedings.pdf|2003]] ([[ftp://gcc.gnu.org/pub/gcc/summit/2003|Individual papers]])
 * [[attachment:2004-GCC-Summit-Proceedings.pdf|2004]] ([[ftp://gcc.gnu.org/pub/gcc/summit/2004|Individual papers]])
 * [[attachment:2005-GCC-Summit-Proceedings.pdf|2005]]
 * [[attachment:2006-GCC-Summit-Proceedings.pdf|2006]]
 * [[attachment:GCC2007-Proceedings.pdf|2007]]
 * [[http://ols.fedoraproject.org/GCC/Reprints-2008/GCC-2008-Proceedings.pdf|2008]] ([[http://ols.fedoraproject.org/GCC/Reprints-2008/|Individual papers]])
 * [[attachment:2009-GCC-Summit-Proceedings.pdf|2009]]
 * [[summit2010|2010]]

== GCC as a research compiler ==
 * [[GCC_Research|Notes for beginners]]

== GCC Research Opportunities Workshop Proceedings ==
 * [[GREPS-2007]] ([[http://sysrun.haifa.il.ibm.com/hrl/greps2007/|GREPS'07 website]])
 * [[GROW-2009]] ([[http://www.doc.ic.ac.uk/~phjk/GROW09/|GROW'09 website]])
 * [[GROW-2010]] ([[http://ctuning.org/workshop-grow10|GROW'10 website]])

== Miscellaneous ==
 * [[History|History of GCC]]
 * [[People]]
 * [[CppConventions|Proposed C++ coding conventions]]
 * [[VolatileAccessComparison|Compiler comparison regarding volatile accesses]]
 * [[Planet_GCC|GCC developer blogs]]
 * [[GCC_glossary|GCC glossary]] and [[abbreviations_and_acronyms|abbreviations and acronyms]]
 * [[DeadlySins|Deadly sins for a compiler writer]]
 * [[Links]]
 * ListOfCompilerBooks
 * OrphanedPages: these pages should be linked from another page, or integrated into other pages and deleted. This doesn't apply to automatic redirections or homepages.
                BESSI1=X*(P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7))))))
        ELSE
                AX=DABS(X)
                Y=3.75/AX
                BESSI1=(DEXP(AX)/DSQRT(AX))*(Q1+Y*(Q2+Y*(Q3+Y*(Q4
     * +Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9))))))))
                IF(X.LT.0.)BESSI1=-BESSI1
        ENDIF
        RETURN
        END
c
C-------------------------------------------------------------------------
 FUNCTION BESSI(N,X)
 IMPLICIT REAL*8 (A-H,O-Z)
C
C FUNCTION RETURNS THE MODIFIED BESSEL FUNCTION OF ORDER N
C FOR ANY REAL X AND N >= 2
C
 IACC=40
 BIGNO=1.0E10
 BIGNI=1.0E-10
 IF (N.LT.2) PAUSE ' bad argument N in BESSI'
 IF (X.EQ.0.0) THEN
  BESSI=0.0
 ELSE
  TOX=2.0/DABS(X)
  BIP=0.0
  BI=1.0
  BESSI=0.
  M=2.*(N+DINT(DSQRT(DBLE(IACC*N))))
  DO 11 J=M,1,-1
   BIM=BIP+DBLE(J)*TOX*BI
   BIP=BI
   BI=BIM
   IF(DABS(BI).GT.BIGNO) THEN
      BESSI=BESSI*BIGNI
      BI=BI*BIGNI
      BIP=BIP*BIGNI
   ENDIF
   IF(J.EQ.N) BESSI=BIP
11 CONTINUE
  BESSI=BESSI*BESSI0(X)/BI
  IF (X.LT.0.0.AND.MOD(N,2).EQ.1) BESSI=-BESSI
 ENDIF
 RETURN
 END


C 19-TH February, 1992 C C INPUT R* AND V INSTEAD MU AND BETA C TRY TO FIND THE CRITICAL GEOMETRIC RATIO FOR RUBBER C C EIGENVALUES SURFACE OF ANTISYMMETRIC DIFFUSE MODES C ************ FOR THE EI SUBREGIME *************** C C FOR INCOMPRESSIBLE ISOTROPIC PRESSURE-INSENSITIVE MATERIAL (RUBBER) C UNDER NO CONFINING STRESSES C C TRY TO COMPARE WITH BEATTY'S (1976) EXPERIMENTAL RESULTS C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C C FOR PLOTTING ALL THOSE ELLIPTIC EIGENVALUES MODE C WITH CHANGING SZZ/(2.*GL) RATIO AND GAMMA (WAVELENGTH) C C DEVELOPED BY C K. T. CHAU C C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C

  • PROGRAM QAMPAPER IMPLICIT REAL*8 (A-H,O-Z) DIMENSION XB1(10),XB2(10) EXTERNAL FEIS,FEIA,ZBRENT CALL CTIME(0)

C C HAVE TO INPUT FIVE MATERIAL PARAMETERS C AND ONE PARAMETER IS CHANGING WITH CURVES C C READ(5,10001)RG,RGK,PR,RR,SRE C READ(5,10002)REG C10001 FORMAT(5F10.0) C10002 FORMAT(F10.0) C C WE HAVE

  • RG =1.0 RGK=0.0 PR =0.5 RR =1.0 SRE=0.0 REG=3.0

C IN THIS PROGRAM C C DEFINE THE FOLLOWING RATIO : C C PR = POISSON RATIO C R = PRESSURE-SENSITIVE COEFFICIENT C RG = GT/GL C RGK = GL/K C SRE = SRR/E C SZE = SZZ/E (THIS THE UNKNOWN TO BE SOLVED) C RVV = R/2V (NON-NORMALITY PARAMETER) C WHERE C GL = LONGITUDINAL SHEAR MODULUS C GT = TRANSVERSE SHEAR MODULUS C E = INSTANTANEOUS TANGENT MODULUS C K = IN-PLANE BULK MODULUS C RR = PRESSURE-SENSITIVITY FACTOR C PR = EF C EFFCTIVE POISSON'S RATIO (V) C

  • DO 100 KK1=1,2

C C KK1=1 FOR SYMMETRIC MODE C KK2=2 FOR ANTISYMMETRIC MODE C

  • WRITE(6,20008)RG,RGK,PR,RR,REG,SRE

C 20008 FORMAT(5X,' APPLIED COEFFICENT OF THE PROBLEM'/

  • 5X,' GT/GL (IN-PLANE BULK MODULUS) =',E12.3/
  • 5X,' GL/K (INSTANTANEOUS E) =',E12.3/
  • 5X,' PR (POISSON RATIO ) =',E12.3/
  • 5X,' RR (PRESSURE-SENSITIVITY) =',E12.3/
  • 5X,' E/GL (BETA INRR MODEL) =',E12.3/
  • 5X,' SRR/E (RADIAL STRESS ) =',E12.3//)

C C

  • TOL=1.E-6

C C GAMMA = N(PI)A/L C WHERE C N = WAVE NUMBER C A = RADIUS OF THE CYLINDER C L = LENGTH OF THE CYLINDER C (GAMMA WILL BE THE X-AXIS OF THE PLOT) C

  • GAMMA=50.0 DG=GAMMA/100.

C C

  • WRITE(6,2011)REG

2011 FORMAT(1X,' E/GL =',E20.6//) C

  • DO 200 KK2=1,1

C C KK2 =1 FOR S<0 C KK2 =2 FOR S>0 C C

  • WRITE(6,292) WRITE(0,292)

292 FORMAT(1X,' GAM SZZ/2GL',//)

  • DO 221 IIB=1,101 GAM=(IIB-1)*DG

C C ELLIPTIC COMPLEX Regime C

  • X1=0.00
    • X2=-0.99 IF(KK2.EQ.2) X2=0.99 N=10

C C NB = NO. OF ROOTS WANTED C N= NO. DIVISION INPUT C

  • NB=2
  • IF(KK1.EQ.1) THEN
    • CALL ZBRAK(FEIS,X1,X2,GAM,XB1,XB2,N,NB)
    ELSE
    • CALL ZBRAK(FEIA,X1,X2,GAM,XB1,XB2,N,NB)
    ENDIF

C C

  • IF(NB.EQ.0) THEN
  • WRITE(0,12011)GAM

12011 FORMAT(1X,'THERE IS NO ROOT (EC REGIME) AT GAM=',F20.5)

  • GO TO 221
  • ENDIF

C XACC=1.E-6

  • DO 776 LL1=1,NB
    • XX1=XB1(LL1) XX2=XB2(LL1)
    IF(KK1.EQ.1) THEN
    • SZZ2=ZBRENT(FEIS,XX1,XX2,GAM)
    ELSE
    • SZZ2=ZBRENT(FEIA,XX1,XX2,GAM)
    ENDIF

C

  • WRITE(6,21111)GAM,SZZ2 WRITE(0,21111)GAM,SZZ2

776 CONTINUE 21111 FORMAT(2E20.6) C 221 CONTINUE C C C Elliptic-Imaginary C C C


C END OF THE EGL (VARYING) LOOP C


200 CONTINUE 100 CONTINUE 999 CONTINUE

  • CALL CTIME (-1) STOP END

C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C DEVELOPED BY C K. T. CHAU C REF: NUMERICAL RECIPES,PP.252 C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C Using Brent's method, find the root of a function FUNC knowm to c lie between X1 and X2. The root, returned as ZBREET, will be c refined until its accuracy is TOL. C

  • FUNCTION ZBRENT(FUNCC,X3,X4,GAM) IMPLICIT REAL*8(A-H,O-Z) EXTERNAL FUNCC ITMAX=250 TOL=1.E-6 EPS=1.E-16 A=X3 B=X4

C

  • FA=FUNCC(A,GAM) FB=FUNCC(B,GAM)

CC C

  • IF (FB*FA.GT.0.) PAUSE 'ROOT MUST BE BRACKETED FOR ZBRENT!' FC=FB DO 11 ITER=1,ITMAX
    • IF (FB*FC.GT.0.0) THEN

C C RENAME A, B AND C AND ADJUST BOUNDING INTERVAL D. C

  • C=A FC=FA D=B-A E=D ENDIF IF (DABS(FC).LT.DABS(FB)) THEN A=B B=C C=A FA=FB FB=FC FC=FA ENDIF

CC C CONVERGENCE CHECK C

  • TOL1=2.*EPS*DABS(B)+0.5*TOL XM=0.5*(C-B) IF (DABS(XM).LE.TOL1.OR.FB.EQ.0.0) THEN ZBRENT=B RETURN ENDIF IF (DABS(E).GE.TOL1.AND.DABS(FA).GT.DABS(FB)) THEN

CC ATTEMPT INVERSE QUADRATIC INTERPOLATION C

  • S=FB/FA IF (A.EQ.C) THEN P=2.*XM*S Q=1.-S ELSE Q=FA/FC R=FB/FC P=S*(2.*XM*Q*(Q-R)-(B-A)*(R-1.)) Q=(Q-1.)*(R-1.)*(S-1.) ENDIF

C C CHECK WHETHER IN BOUNDS C

  • IF (P.GT.0.0) Q=-Q P=DABS(P) IF (2.*P.LT.DMIN1(3.*XM*Q-DABS(TOL1*Q),DABS(E*Q))) THEN

C C ACCEPT INTERPOLATION C

  • E=D D=P/Q ELSE

C C INTERPOLATION FAILED, USE BISECTION. C

  • D=XM E=D ENDIF ELSE

C C BOUNDS DECREASES TOO SLOWLY, USE BISECTION. C

  • D=XM E=D ENDIF

C C MOVE LAST BEST GUESS TO A C

  • A=B FA=FB IF(DABS(D).GT.TOL1) THEN

C C EVALUATE NEW TRIAL ROOT C

  • B=B+D ELSE B=B+DSIGN(TOL1,XM) ENDIF

C

  • FB=FUNCC(B,GAM)

C 11 CONTINUE

  • PAUSE 'ZBRENT EXCEEDING MAXIMUM ITERATIONS' ZBRENT=B RETURN END

C C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C SUBROUTINE TO BRACKET THE INTERVAL FROM X1 TO X2 C

  • SUBROUTINE ZBRAK(FX,X1,X2,GAM,XBB1,XBB2,N,NB)

C C GIVEN A FUNCTION FX DEFINED ON THE INTERVAL FROM X1 TO X2 C SUBDIVIDE THE INTERVAL INTO N EQUALLY SPACED SEGMENTS, AND C SEARCH FOR ZERO CROSSINGS OF THE FUNCTION. NB IS INPUT AS C THE MAXIMUM NUMBER OF ROOTS SOUGHT, AND IS RESET TO THE C NUMBER OF BRACKETING PAIRS XB1, XB2 THAT ARE FOUND. C C

  • IMPLICIT REAL*8 (A-H,O-Z) DIMENSION XBB1(10),XBB2(10) EXTERNAL FX NBB=NB NB=0 X=X1 DX=(X2-X1)/N

C

  • FP=FX(X,GAM)

C C LOOP OVER ALL INTERVALS C

  • DO 11 I=1,N
    • X=X+DX

C

  • FC=FX(X,GAM)

C

  • IF(FC*FP.LT.0.0) THEN

C IF SIGN CHANGE OCCURS THEN RECORD VALUES FOR THE BOUNDS C

  • NB=NB+1 XBB1(NB)=X-DX XBB2(NB)=X ENDIF FP=FC IF (NBB.EQ.NB) RETURN

11 CONTINUE

  • RETURN END

C---|


1


|


2


|


3


|


4


|


5


|


6


|


712 C

  • FUNCTION FEIA(SZZ,GAM) IMPLICIT REAL*8 (A-H,O-Z) EXTERNAL BESSI0,BESSI1,BESSI

CC C NEW SUBROUTINE FOR ANTI-SYMMETRIC MODES C

  • SS=SZZ VA=(1.-SS) VB=2. VC=(1.+SS)

C V3 FOR THE SECOND DIFFERENTIAL EQUATION

  • U3=DSQRT(1.+SS) U32=U3*U3

C C Elliptic-Complex Regime C C IF((VB*VB-4.*VA*VC).LT.0.0.AND.(VA*VC).GT.0.0) THEN C IF(VA.EQ.0.0) PAUSE 'A=0 IN FECA' C C CHECK POSITIVENESS OF P AND Q C C CALCULATION OF P AND Q

  • PP=DSQRT((VB+DSQRT(VB*VB-4.*VA*VC))/(2.*VA)) QQ=DSQRT((VB-DSQRT(VB*VB-4.*VA*VC))/(2.*VA)) PP2=PP*PP QQ2=QQ*QQ GP=GAM*PP GQ=GAM*QQ

C

  • BIP0=BESSI0(GP) BIP1=BESSI1(GP) BIP2=BESSI(2,GP) BIP3=BESSI(3,GP) BIQ0=BESSI0(GQ) BIQ1=BESSI1(GQ) BIQ2=BESSI(2,GQ) BIQ3=BESSI(3,GQ)

C

  • ARG=GAM*U3 BI0 =BESSI0(ARG) BI1 =BESSI1(ARG) BI2 =BESSI(2,ARG) BI3 =BESSI(3,ARG)

C---|


1


|


2


|


3


|


4


|


5


|


6


|


712

  • Y1 = -(1.+SS) Y2 = Y1 Y3 = 2.0 Y4 = -(1.-SS) Y5 = (1.-SS) Y6 = Y5 Y7 = 1.0 Y8 = Y7

C

  • IF(GAM.EQ.0.0) GO TO 777 WW1 = -4.*U32*Y3*BI2/ARG WW2 = U3*Y6*BI1/ARG WW3 = U32/4.*(2.*Y7/ARG*(BI0+BI2-2.*BI1/ARG)-Y8*(3.*BI1+BI3))

C

  • VPP=((4.*Y1+3.*Y3)*PP2-4.*Y2)*BIP1+PP2*Y3*BIP3 VQP=0.5*PP*(Y4*PP2-Y5)*(BIP0+BIP2) VRP=PP2*Y3*BIP2/GP VPQ=((4.*Y1+3.*Y3)*QQ2-4.*Y2)*BIQ1+QQ2*Y3*BIQ3 VQQ=0.5*QQ*(Y4*QQ2-Y5)*(BIQ0+BIQ2) VRQ=QQ2*Y3*BIQ2/GQ

C C---|


1


|


2


|


3


|


4


|


5


|


6


|


712

  • GO TO 778

777 FEIA=U32*Y3*Y3*Y5+4.*Y2*Y3*Y6+U32*(Y7-3.*Y8)*(Y5*(4.*Y1+3.*Y3)

  • -4.*Y2*Y4)
    • GO TO 999

C ELSE 778 FEIA=WW1*(VQP*VRQ-VQQ*VRP)+WW2*(VPQ*VRP-VPP*VRQ)

  • +WW3*(VPP*VQQ-VPQ*VQP)

999 RETURN

  • END

C


C EIGENVALUE EQUATION FOR SYMMETRIC GEOMETRIC DIFFUSE MODE C


  • FUNCTION FEIS(SZZ,GAM) IMPLICIT REAL*8 (A-H,O-Z) EXTERNAL BESSI1,BESSI0,BESSI

C

  • SS=SZZ VA=(1.-SS) VB=2. VC=(1.+SS)

C C Elliptic-Imaginary Subregime C

  • PP=DSQRT((VB+DSQRT(VB*VB-4.*VA*VC))/(2.*VA)) QQ=DSQRT((VB-DSQRT(VB*VB-4.*VA*VC))/(2.*VA))

c PP=DSQRT((2.*DSQRT(VA*VC)-VB)/(4.*VA)) c QQ=DSQRT((2.*DSQRT(VA*VC)+VB)/(4.*VA))

  • PP2=PP*PP QQ2=QQ*QQ GP=GAM*PP GQ=GAM*QQ

C

  • BIP0=BESSI0(GP) BIP1=BESSI1(GP) BIP2=BESSI(2,GP) BIQ0=BESSI0(GQ) BIQ1=BESSI1(GQ) BIQ2=BESSI(2,GQ)

C---|


1


|


2


|


3


|


4


|


5


|


6


|


712

  • Y1 = -(1.+SS) Y2 = Y1 Y3 = 2.0 Y4 = -(1.-SS) Y5 = (1.-SS) Y6 = Y5 Y7 = 1.0 Y8 = Y7

C

  • VPP=2.*((2.*Y1+Y3)*PP2-2.*Y2)*BIP0+PP2*Y3*2.*BIP2 VPQ=2.*((2.*Y1+Y3)*QQ2-2.*Y2)*BIQ0+QQ2*Y3*2.*BIQ2 VQP=2.*PP*(Y4*PP2-Y5)*BIP1 VQQ=2.*QQ*(Y4*QQ2-Y5)*BIQ1

C C

  • IF(GAM.EQ.0.0) THEN FEIS=VC*Y4*(2.*Y1+Y3)+2.*Y2*(VA*Y5-VB*Y4) GO TO 999 ELSE FEIS=VPP*VQQ-VPQ*VQP ENDIF

999 RETURN

  • END

C C SUBROUTINE FOR TIME CALCULATION FOR THE RUNNING TIME C **************************************************** C

  • INTERFACE TO SUBROUTINE TIME (N,STR) CHARACTER*10 STR [NEAR,REFERENCE] INTEGER*2 N [VALUE] END

$LARGE

  • SUBROUTINE CTIME(N) IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*1 C1,C2,C4,C5,C7,C8 CHARACTER*10 CLOCK COMMON/CPUT/CT,PT CALL TIME (10,CLOCK) C1=CLOCK(1:1) C2=CLOCK(2:2) C4=CLOCK(4:4) C5=CLOCK(5:5) C7=CLOCK(7:7) C8=CLOCK(8:8) HOUR=(ICHAR(C1)-48)*10+ICHAR(C2)-48 MINU=(ICHAR(C4)-48)*10+ICHAR(C5)-48 SECO=(ICHAR(C7)-48)*10+ICHAR(C8)-48 A=HOUR*3600.+MINU*60.+SECO IF ( N.EQ.0) GO TO 100 IF ( N.LT.0) GO TO 200 IF ( N.GT.0) GO TO 300

300 CT = A

  • IF (CT.LT.PT) CT = CT + 3600.*24. T = CT - PT WRITE (6,1000) T WRITE (0,1000) T

1000 FORMAT (41X, ' CURRENT TIME (SEC.)>>>>' ,F10.0,/)

  • RETURN

100 PT = A

  • ST = 0.0 WRITE (6,2000) ST WRITE (0,2000) ST

2000 FORMAT (41X, ' STARTING TIME (SEC.)>>>>',F10.0,/)

  • RETURN

200 CT = A

  • IF (CT.LT.PT) CT = CT + 3600.*24. T = CT - PT WRITE (6,3000) T WRITE (0,3000) T

3000 FORMAT (41X, ' ENDING TIME (SEC.)>>>>',F10.0,/)

  • RETURN END

C---|


1


|


2


|


3


|


4


|


5


|


6


|


712

  • FUNCTION BESSI0(X) IMPLICIT REAL*8 (A-H,O-Z) DATA P1,P2,P3,P4,P5,P6,P7/1.D0,3.5156229D0,3.0899424D0,
  • 1.2067492D0,0.2659732D0,0.360768D-1,0.45813D-2/
    • DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9/0.39894228D0,0.1328592D-1,
  • 0.225319D-2,-0.157565D-2,0.916281D-2,-0.2057706D-1,
  • 0.2635537D-1,-0.1647633D-1,0.392377D-2/
    • IF(DABS(X).LT.3.75) THEN
      • Y=(X/3.75)**2. BESSI0=P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7)))))
      ELSE
      • AX=DABS(X) Y=3.75/AX BESSI0=(DEXP(AX)/DSQRT(AX))*(Q1+Y*(Q2+Y*(Q3+Y*(Q4
  • +Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9))))))))
    • ENDIF RETURN END

C


  • FUNCTION BESSI1(X) IMPLICIT REAL*8 (A-H,O-Z)

C---|


1


|


2


|


3


|


4


|


5


|


6


|


712

  • DATA P1,P2,P3,P4,P5,P6,P7/0.5D0,0.87890594D0,0.51498869D0,
  • 0.15084934D0,0.2658733D-1,0.301532D-2,0.32411D-3/
    • DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9/0.39894228D0,-0.3988024D-1,
  • -0.362018D-2,0.163801D-2,-0.1031555D-1,0.2282967D-1,
  • -0.2895312D-1,0.1787654D-1,-0.420059D-2/
    • IF(DABS(X).LT.3.75) THEN
      • Y=(X/3.75)**2. BESSI1=X*(P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7))))))
      ELSE
      • AX=DABS(X) Y=3.75/AX BESSI1=(DEXP(AX)/DSQRT(AX))*(Q1+Y*(Q2+Y*(Q3+Y*(Q4
  • +Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9))))))))
    • IF(X.LT.0.)BESSI1=-BESSI1
    • ENDIF RETURN END

c C


  • FUNCTION BESSI(N,X) IMPLICIT REAL*8 (A-H,O-Z)

C C FUNCTION RETURNS THE MODIFIED BESSEL FUNCTION OF ORDER N C FOR ANY REAL X AND N >= 2 C

  • IACC=40 BIGNO=1.0E10 BIGNI=1.0E-10 IF (N.LT.2) PAUSE ' bad argument N in BESSI' IF (X.EQ.0.0) THEN
    • BESSI=0.0
    ELSE
    • TOX=2.0/DABS(X) BIP=0.0 BI=1.0 BESSI=0. M=2.*(N+DINT(DSQRT(DBLE(IACC*N)))) DO 11 J=M,1,-1
      • BIM=BIP+DBLE(J)*TOX*BI BIP=BI BI=BIM IF(DABS(BI).GT.BIGNO) THEN
        • BESSI=BESSI*BIGNI BI=BI*BIGNI BIP=BIP*BIGNI
        ENDIF IF(J.EQ.N) BESSI=BIP

11 CONTINUE

  • BESSI=BESSI*BESSI0(X)/BI IF (X.LT.0.0.AND.MOD(N,2).EQ.1) BESSI=-BESSI
  • ENDIF RETURN END



None: HomePage (last edited 2018-09-21 13:49:53 by IainSandoe)