diff --git a/wafo/source/mreg/intfcmod.f b/wafo/source/mreg/intfcmod.f index 9cda573..927dd3e 100644 --- a/wafo/source/mreg/intfcmod.f +++ b/wafo/source/mreg/intfcmod.f @@ -145,8 +145,8 @@ C Covariance matrices COV1=r'(T-T), COV2=r''(T-T) and COV3=r'''(T-T) C Dimension of COV1, COV2 should be atleast N*N. C USE SIZEMOD -! IMPLICIT NONE -C INTEGER, PARAMETER:: NMAX = 101, RDIM = 10201 +! IMPLICIT NONE +C INTEGER, PARAMETER:: NMAX = 101, RDIM = 10201 REAL*8, PARAMETER:: ZERO = 0.0d0 REAL*8, intent(inout) :: XL0,XL2,XL4 REAL*8, DIMENSION(N,5), intent(in) :: COV diff --git a/wafo/source/mvnprd/mvnprd.f b/wafo/source/mvnprd/mvnprd.f index 03f97cc..be1de64 100644 --- a/wafo/source/mvnprd/mvnprd.f +++ b/wafo/source/mvnprd/mvnprd.f @@ -1,12 +1,12 @@ -C +C C f2py -m mvnprd -h mvnprd.pyf mvnprd.f only: mvnprd C edit mvnprd.pyf with input and output and then C f2py mvnprd.pyf mvnprd.f -c --fcompiler=gnu95 --compiler=mingw32 -lmsvcr71 C -C f2py --fcompiler=gnu95 --compiler=mingw32 -lmsvcr71 -m mvnprd -c mvnprd.f +C f2py --fcompiler=gnu95 --compiler=mingw32 -lmsvcr71 -m mvnprd -c mvnprd.f C Altarnative: compile mvnprd and link to it through mvnprd_interface.f -C +C C gfortran -fPIC -c mvnprd.f C f2py -m mvnprdmod -c mvnprd.o mvnprd_interface.f --fcompiler=gnu95 --compiler=mingw32 -lmsvcr71 C @@ -18,7 +18,7 @@ C f2py -m mvnprdmod -c mvnprd.obj mvnprd_interface.f --fcompiler=compaqv --comp ! ,MVNPRD and MVSTUD subroutines for computing multivariate normal ! or student T probabilities with product correlation structure. The ! file should compile without errors on (Fortran77) standard Fortran -! compilers. +! compilers. * * The mex-interface was written by * Per Andreas Brodtkorb @@ -45,7 +45,7 @@ C Tel.: (905) 525-9140 (Ext. 27104) * RHO REAL, array of coefficients defining the correlation * coefficient by: * correlation(I,J) = RHO(I)*RHO(J) for J/=I -* where +* where * 1 < RHO(I) < 1 * A REAL, array of lower integration limits. * B REAL, array of upper integration limits. @@ -66,13 +66,13 @@ C Tel.: (905) 525-9140 (Ext. 27104) * INFORM INTEGER, termination status parameter: * 0, if normal completion with ERROR < EPS; * 1, if N > 100 or N < 1. -* 2, IF any abs(rho)>=1 +* 2, IF any abs(rho)>=1 * 4, if ANY(B(I)<=A(i)) * 5, if number of terms computed exceeds maximum number of * evaluation points * 6, if fault accurs in normal subroutines * 7, if subintervals are too narrow or too many -* 8, if bounds exceeds abseps +* 8, if bounds exceeds abseps * * * MVNPRDMEX calculates multivariate normal or student T probability @@ -836,7 +836,7 @@ C DOUBLE PRECISION :: HNC,PROB,BND INTEGER :: NN, MAXDF,I,IERC,NDF,N,IFLT PARAMETER (NN=100, MAXDF = 150) - integer :: INF(*) + integer :: INF(*) DOUBLE PRECISION :: A(*),B(*),BPD(*),D(*),F(3), & AA(NN),BB(NN) DOUBLE PRECISION :: ERB2, ERRB, AX,BX,XX @@ -855,7 +855,7 @@ C ENDIF BND = ZERO IFLT = 0 - + ERB2 = ERRB C C CHECK IF D.F. EXCEED MAXDF; IF YES, THEN PROB @@ -1086,7 +1086,7 @@ C RETURN END -C * * * * * * * * * * * * * * * * * * * * * * * * * * * * +C * * * * * * * * * * * * * * * * * * * * * * * * * * * * C Charles Dunnett C Dept. of Mathematics and Statistics C McMaster University diff --git a/wafo/source/mvnprd/mvnprodcorrprb.f b/wafo/source/mvnprd/mvnprodcorrprb.f index ae000cc..02b621a 100644 --- a/wafo/source/mvnprd/mvnprodcorrprb.f +++ b/wafo/source/mvnprd/mvnprodcorrprb.f @@ -304,11 +304,11 @@ C DOUBLE PRECISION, DIMENSION(6) :: P C------------------------------------------------------------------ C Coefficients for approximation to erf in first interval C------------------------------------------------------------------ - DOUBLE PRECISION, PARAMETER, DIMENSION(5) :: + DOUBLE PRECISION, PARAMETER, DIMENSION(5) :: & A = (/ 3.16112374387056560D00, & 1.13864154151050156D02,3.77485237685302021D02, & 3.20937758913846947D03, 1.85777706184603153D-1/) - DOUBLE PRECISION, PARAMETER, DIMENSION(4) :: + DOUBLE PRECISION, PARAMETER, DIMENSION(4) :: & B = (/2.36012909523441209D01,2.44024637934444173D02, & 1.28261652607737228D03,2.84423683343917062D03/) C------------------------------------------------------------------ @@ -320,7 +320,7 @@ C------------------------------------------------------------------ 2 8.81952221241769090D02,1.71204761263407058D03, 3 2.05107837782607147D03,1.23033935479799725D03, 4 2.15311535474403846D-8/) - DOUBLE PRECISION, DIMENSION(8) :: + DOUBLE PRECISION, DIMENSION(8) :: & D =(/1.57449261107098347D01,1.17693950891312499D02, 1 5.37181101862009858D02,1.62138957456669019D03, 2 3.29079923573345963D03,4.36261909014324716D03, @@ -328,14 +328,14 @@ C------------------------------------------------------------------ C------------------------------------------------------------------ C Coefficients for approximation to erfc in third interval C------------------------------------------------------------------ - DOUBLE PRECISION, parameter, + DOUBLE PRECISION, parameter, & DIMENSION(6) :: P =(/3.05326634961232344D-1, & 3.60344899949804439D-1, 1 1.25781726111229246D-1,1.60837851487422766D-2, 2 6.58749161529837803D-4,1.63153871373020978D-2/) - DOUBLE PRECISION, parameter, + DOUBLE PRECISION, parameter, & DIMENSION(5) :: Q =(/2.56852019228982242D00, - & 1.87295284992346047D00, + & 1.87295284992346047D00, 1 5.27905102951428412D-1,6.05183413124413191D-2, 2 2.33520497626869185D-3/) C------------------------------------------------------------------ @@ -3090,29 +3090,29 @@ c ! Subdivide the interval and create two new vectors in the stack, ! one of which overwrites the vector just processed. ! -! v(:,k) = [fx1,fx2,fx3,fx4,fx5,x1,h,S,SL,SR] - kp1 = k + 1; -! Process right interval - v(1,kp1) = v(3,k); !fx1R - v(2,kp1) = fx(3); !fx2R - v(3,kp1) = v(4,k); !fx3R - v(4,kp1) = fx(4); !fx4R - v(5,kp1) = v(5,k); !fx5R - v(6,kp1) = v(6,k) + four * h; ! x1R - v(7,kp1) = h; - v(8,kp1) = v(10,k); ! S - v(9:10,kp1) = Sn(3:4); ! SL, SR -! Process left interval - v(5,k) = v(3,k); ! fx5L - v(4,k) = fx(2); ! fx4L - v(3,k) = v(2,k); ! fx3L - v(2,k) = fx(1); ! fx2L -! v(1,k) unchanged fx1L -! v(6,k) unchanged x1L - v(7,k) = h; - v(8,k) = v(9,k); ! S - v(9:10,k) = Sn(1:2); ! SL, SR - k = kp1; +! v(:,k) = [fx1,fx2,fx3,fx4,fx5,x1,h,S,SL,SR] + kp1 = k + 1; +! Process right interval + v(1,kp1) = v(3,k); !fx1R + v(2,kp1) = fx(3); !fx2R + v(3,kp1) = v(4,k); !fx3R + v(4,kp1) = fx(4); !fx4R + v(5,kp1) = v(5,k); !fx5R + v(6,kp1) = v(6,k) + four * h; ! x1R + v(7,kp1) = h; + v(8,kp1) = v(10,k); ! S + v(9:10,kp1) = Sn(3:4); ! SL, SR +! Process left interval + v(5,k) = v(3,k); ! fx5L + v(4,k) = fx(2); ! fx4L + v(3,k) = v(2,k); ! fx3L + v(2,k) = fx(1); ! fx2L +! v(1,k) unchanged fx1L +! v(6,k) unchanged x1L + v(7,k) = h; + v(8,k) = v(9,k); ! S + v(9:10,k) = Sn(1:2); ! SL, SR + k = kp1; endif enddo ! while if (epsi0) iflg = ior(iflg,kflg) end do - if (epsi