git-svn-id: https://lxsd.femto-st.fr/svn/fvn@33 b657c933-2333-4658-acf2-d3c7c2708721

daniau
1 parent c90ee9d32e
Showing 1 changed file with 14 additions and 14 deletions Side-by-side Diff
fvnlib.f90
@@ -2127,18 +2127,18 @@
 !
 subroutine fvn_z_muller (f,eps,eps1,kn,nguess,n,x,itmax,infer,ier)
      implicit none
-      double precision :: rzero,rten,rhun,rp01,ax,eps1,qz,eps,tpq
+      double precision :: rzero,rten,rhun,rp01,ax,eps1,qz,eps,tpq,eps1w
       double complex ::   d,dd,den,fprt,frt,h,rt,t1,t2,t3, &
                           tem,z0,z1,z2,bi,xx,xl,y0,y1,y2,x0, &
                           zero,p1,one,four,p5
-      
+
       double complex, external :: f
       integer :: ickmax,kn,nguess,n,itmax,ier,knp1,knpn,i,l,ic, &
                     knpng,jk,ick,nn,lm1,errcode
       double complex :: x(kn+n)
       integer :: infer(kn+n)
-      
-      
+
+
       data                zero/(0.0,0.0)/,p1/(0.1,0.0)/, &
                           one/(1.0,0.0)/,four/(4.0,0.0)/, &
                           p5/(0.5,0.0)/, &
@@ -2150,8 +2150,8 @@
                 return
             end if
             !eps1 = rten **(-nsig)
-            eps1 = min(eps1,rp01)
-            
+            eps1w = min(eps1,rp01)
+
             knp1 = kn+1
             knpn = kn+n
             knpng = kn+nguess
@@ -2199,7 +2199,7 @@
                         exit icloop
                     end if
  
-   
+ 
                     z2 = fprt
                     y2 = frt
 !                                  begin main algorithm
@@ -2241,7 +2241,7 @@
                         h = d*h
                         rt = rt + h
 !                                  check convergence of the first kind  
-                        if (cdabs(h) .le. eps1*max(cdabs(rt),ax)) then
+                        if (cdabs(h) .le. eps1w*max(cdabs(rt),ax)) then
                             if (ic .ne. 0) then
                                 exit icloop
                             end if
  
  
  
@@ -2289,21 +2289,21 @@
         infer(l) = jk
         l = l+1
       end do rloop
-      
+
       contains
         subroutine trans_rt()
-           tem = rten*eps1
+           tem = rten*eps1w
            if (cdabs(rt) .gt. ax) tem = tem*rt
            rt = rt+tem
            d = (h+tem)*d/h
            h = h+tem
         end subroutine trans_rt
-        
+ 
         subroutine deflated_work(errcode)
             ! errcode=0 => no errors
             ! errcode=1 => jk>itmax or convergence of second kind achieved
             integer :: errcode,flag
-            
+ 
             flag=1
     loop1:  do while(flag==1)
                 errcode=0
  
@@ -2334,9 +2334,9 @@
                 errcode=1
                 return
             end if
-            
+
         end subroutine deflated_work
-      
+
       end subroutine
...	...	@@ -2127,18 +2127,18 @@
2127	2127	!
2128	2128	subroutine fvn_z_muller (f,eps,eps1,kn,nguess,n,x,itmax,infer,ier)
2129	2129	implicit none
2130		- double precision :: rzero,rten,rhun,rp01,ax,eps1,qz,eps,tpq
	2130	+ double precision :: rzero,rten,rhun,rp01,ax,eps1,qz,eps,tpq,eps1w
2131	2131	double complex :: d,dd,den,fprt,frt,h,rt,t1,t2,t3, &
2132	2132	tem,z0,z1,z2,bi,xx,xl,y0,y1,y2,x0, &
2133	2133	zero,p1,one,four,p5
2134		-
	2134	+
2135	2135	double complex, external :: f
2136	2136	integer :: ickmax,kn,nguess,n,itmax,ier,knp1,knpn,i,l,ic, &
2137	2137	knpng,jk,ick,nn,lm1,errcode
2138	2138	double complex :: x(kn+n)
2139	2139	integer :: infer(kn+n)
2140		-
2141		-
	2140	+
	2141	+
2142	2142	data zero/(0.0,0.0)/,p1/(0.1,0.0)/, &
2143	2143	one/(1.0,0.0)/,four/(4.0,0.0)/, &
2144	2144	p5/(0.5,0.0)/, &
...	...	@@ -2150,8 +2150,8 @@
2150	2150	return
2151	2151	end if
2152	2152	!eps1 = rten **(-nsig)
2153		- eps1 = min(eps1,rp01)
2154		-
	2153	+ eps1w = min(eps1,rp01)
	2154	+
2155	2155	knp1 = kn+1
2156	2156	knpn = kn+n
2157	2157	knpng = kn+nguess
...	...	@@ -2199,7 +2199,7 @@
2199	2199	exit icloop
2200	2200	end if
2201	2201
2202		-
	2202	+
2203	2203	z2 = fprt
2204	2204	y2 = frt
2205	2205	! begin main algorithm
...	...	@@ -2241,7 +2241,7 @@
2241	2241	h = d*h
2242	2242	rt = rt + h
2243	2243	! check convergence of the first kind
2244		- if (cdabs(h) .le. eps1*max(cdabs(rt),ax)) then
	2244	+ if (cdabs(h) .le. eps1w*max(cdabs(rt),ax)) then
2245	2245	if (ic .ne. 0) then
2246	2246	exit icloop
2247	2247	end if
2248	2248
2249	2249
2250	2250
...	...	@@ -2289,21 +2289,21 @@
2289	2289	infer(l) = jk
2290	2290	l = l+1
2291	2291	end do rloop
2292		-
	2292	+
2293	2293	contains
2294	2294	subroutine trans_rt()
2295		- tem = rten*eps1
	2295	+ tem = rten*eps1w
2296	2296	if (cdabs(rt) .gt. ax) tem = tem*rt
2297	2297	rt = rt+tem
2298	2298	d = (h+tem)*d/h
2299	2299	h = h+tem
2300	2300	end subroutine trans_rt
2301		-
	2301	+
2302	2302	subroutine deflated_work(errcode)
2303	2303	! errcode=0 => no errors
2304	2304	! errcode=1 => jk>itmax or convergence of second kind achieved
2305	2305	integer :: errcode,flag
2306		-
	2306	+
2307	2307	flag=1
2308	2308	loop1: do while(flag==1)
2309	2309	errcode=0
2310	2310
...	...	@@ -2334,9 +2334,9 @@
2334	2334	errcode=1
2335	2335	return
2336	2336	end if
2337		-
	2337	+
2338	2338	end subroutine deflated_work
2339		-
	2339	+
2340	2340	end subroutine
2341	2341
2342	2342