C August 1990 C Dear Colleague: C The following are the files of the programs 1-5 C as they have been used in our book C "Numerical Continuation Methods", C Springer Series in Computational Mathematics 13. C The exclamation mark ! is a comment symbol C which is not understood by all FORTRAN compilers. C If this is the case for you, please replace it C by the following three ASCII symbols: C RETURN C BLANK C Let us remind you of some disclaimers in the preface: C The programs are primarily to be regarded as illustrations. C We hope that the reader will experiment with our illustration C programs and be led to make improvements and adaptations C suited to his particular applications. C We emphasize that our programs should not be regarded C as perfected library programs. C Thank you for your interest! C Best regards, Eugene L. Allgower and Kurt Georg. C Program 1 program cont !continuation method, follows a curve H(u) = 0 !Euler predictor, Newton-correctors !stepsize control by asymptotic estimates !Jacobian is evaluated only at predictor point !stepsize is monitored by two different values: !1. contraction rate in corrector steps !2. distance to the curve !stops at a point x such that x(n1) = 0 !arrays: parameter(n = 10, n1 = n+1) !dimension of the problem dimension b(n1,n) !transpose of Jacobian dimension q(n1,n1) !orth. matrix for QR dec. of b dimension x(n1), u(n1) !current points on the curve dimension t(n1) !tangent vector dimension y(n) !stores values y := H(x) logical succ, newt !parameters: tol = 1.e-4 !tolerance for corrector iteration ctmax = 0.6 !maximal contr. rate in corrector step dmax = .4 !maximal distance to curve hmax = 1. !maximal stepsize hmin = 1.e-5 !minimal stepsize maxjac = 1000 !maximal number of Jacobian evaluations cdmax = 1.e3 !maximum for condition estimate fmax = 2. !maximal factor for acceleration h = .03 !initial stepsize eta = .1 !perturbation to avoid cancellation !when calculating the contraction rate !main program open(1, file='cont.dat') !output file call stpnt(x, n1) !user defined starting point, H(x) = 0 newt = .false. mapct = 0 !counts the calls of the map H jacct = 0 !counts the calls of the Jacobian $H'$ call jacob(b, x, n, n1) !b := transpose of Jacobian at x jacct = jacct + 1 call decomp(b, q, cond, n, n1) !b, q := orthog. decomp. of b if (cond .gt. cdmax) then write(1,*) ' bad cond. estimate in init. point = ', cond write(*,*) ' bad cond. estimate in init. point = ', cond stop endif do 91 k = 1, n1 !tangent t(k) = q(n1, k) 91 continue call setor(or, t, n1) !set orientation 12 continue !begin PC loop if (abs(h).lt.hmin) then write(1,*) ' failure at minimal stepsize' write(*,*) ' failure at minimal stepsize' stop endif if (jacct .gt. maxjac) then write(*,*) ' maximal number of Jacobian eval. exceeded' write(1,*) ' maximal number of Jacobian eval. exceeded' stop endif do 92 k = 1, n1 u(k) = x(k) + h * or * t(k) !predictor step 92 continue fac = 1./ fmax !initialize deceleration factor call jacob(b, u, n, n1) !b := transpose of Jacobian at u jacct = jacct + 1 call decomp(b, q, cond, n, n1) !decompose b if (cond .gt. cdmax) goto 21 iter = 0 !counts the corrector iterations 93 iter = iter + 1 !begin corrector loop call map(u, y, n, n1) mapct = mapct + 1 call newton(q, b, u, y, dist, n, n1) if (dist.gt.dmax) goto 21 fac = max(fac, sqrt(dist/dmax)*fmax) if (iter.ge.2) then contr = dist / (disto + tol*eta) !contraction rate if (contr.gt.ctmax) goto 21 fac = max(fac, sqrt(contr/ctmax)*fmax) endif if (dist.lt.tol) goto 22 !corrector successful disto = dist goto 93 !end corrector loop 21 h = h / fmax !PC not accepted goto 12 22 continue !PC step accepted succ = .false. if (u(n1).ge.1.) newt = .true. !switch to Newton steplength if (newt) then h = - (u(n1) - 1.) / q(n1, n1) if (abs(h).lt.hmin) succ = .true. !solution point found else if (fac.gt.fmax) fac = fmax h = min(abs(h/fac), hmax) !steplength adaptation if (h.gt.hmax) h = hmax endif do 94 k = 1, n1 x(k) = u(k) !new point on curve t(k) = q(n1, k) !new tangent 94 continue if (succ) then !stopping the curve tracing write(1,*) ' success with', mapct,' calls of ''map'' and', * jacct, ' calls of ''jacob''' write(*,*) ' success with', mapct,' calls of ''map'' and', * jacct, ' calls of ''jacob''' write(1,*) write(*,*) write(1,*) ' solution vector:' write(*,*) ' solution vector:' write(1,*) ' ===============' write(*,*) ' ===============' do 95 k = 1, n write(1,*) ' x(', k, ') = ', x(k) write(*,*) ' x(', k, ') = ', x(k) 95 continue stop endif goto 12 end subroutine map(x, y, n, n1) !user defined !input: x output: y = H(x) !H(x) = 0 defines the curve to be traced dimension x(n1), y(n) s = 0.0 do 91 i = 1, n s = s + x(i) 91 continue do 92 i = 1, n y(i) = x(i) - x(n1) * exp(cos(i * s)) 92 continue return end subroutine jacob(b, x, n, n1) !user defined !input: x output: b !evaluates the transpose b of the Jacobian at x dimension b(n1,n), x(n1) s = 0.0 do 91 i = 1, n s = s + x(i) 91 continue do 92 k = 1, n1 do 93 i = 1, n if (k.eq.n1) then b(k, i) = -exp(cos(i * s)) elseif (i.eq.k) then b(k,i)=1.+x(n1)*exp(cos(i*s))*sin(i*s)*i else b(k,i) = x(n1)*exp(cos(i*s))*sin(i*s)*i endif 93 continue 92 continue return end subroutine stpnt(x, n1) !user defined !output: x = starting point on curve dimension x(n1) do 91 k = 1, n1 x(k) = 0.0 91 continue return end subroutine setor(or, t, n1) !user defined !input: t output: or(t) !decides in which direction the curve will be traversed dimension t(n1) if (t(n1).gt.0.) then or = 1.0 else or = -1.0 endif return end subroutine givens(b, q, c1, c2, l1, l2, l3, n, n1) !input: b, q, c1, c2, l1, l2, l3 !output: b, q, c1, c2 !one Givens rotation is performed --- !on rows l1 and l2 of b and q !the rotation maps c1, c2 onto sqrt(c1**2+c2**2), 0 dimension b(n1, n), q(n1, n1) if (abs(c1)+abs(c2) .eq. 0.) return if (abs(c2) .ge. abs(c1)) then sn = sqrt(1. + (c1/c2)**2) * abs(c2) else sn = sqrt(1. + (c2/c1)**2) * abs(c1) endif s1 = c1/sn s2 = c2/sn do 91 k = 1, n1 sv1 = q(l1, k) sv2 = q(l2, k) q(l1, k) = s1 * sv1 + s2 * sv2 q(l2, k) = -s2 * sv1 + s1 * sv2 91 continue do 92 k = l3, n sv1 = b(l1, k) sv2 = b(l2, k) b(l1, k) = s1 * sv1 + s2 * sv2 b(l2, k) = -s2 * sv1 + s1 * sv2 92 continue c1 = sn c2 = 0.0 return end subroutine decomp(b, q, cond, n, n1) !input: b, output: b, q, cond !a QR decomposition for b is stored in q, b --- !by using Givens rotations on b and q = id --- !until b is upper triangular !a very coarse condition estimate cond is provided dimension b(n1, n), q(n1, n1) do 91 k = 1, n1 !start with q := id do 92 l = 1, n1 q(k, l) = 0.0 92 continue q(k, k) = 1.0 91 continue do 93 m = 1, n !successive Givens transformations do 94 k = m+1, n1 call givens(b, q, b(m, m), b(k, m), m, k, m+1, n, n1) 94 continue 93 continue cond = 0. !very coarse condition estimate do 95 i = 2, n do 96 k = 1, i - 1 cond = max(cond, abs(b(k,i)/b(i,i))) 96 continue 95 continue return end subroutine newton(q, b, u, y, d, n, n1) !input q, b, u, y = H(u), n, n1 !output u, d !y is changed !a Newton step $u := u - A^+ H(u)$ is performed --- !where A approximates the current Jacobian $H'$ !q, b = QR decomposition of $A^*$ !d = length of Newton step dimension q(n1, n1), b(n1, n), u(n1), y(n) do 91 k = 1, n do 92 l = 1, k-1 y(k) = y(k) - b(l, k) * y(l) 92 continue y(k) = y(k) / b(k, k) 91 continue d = 0. do 93 k = 1, n1 s = 0.0 do 94 l = 1, n s = s + q(l, k) * y(l) 94 continue u(k) = u(k) - s d = d + s**2 93 continue d = sqrt(d) return end C Program 2 program plhom !Piecewise linear homotopy method !in the sense of Eaves and Saigal, see section 13.4 !with automatic pivoting steps and tentative Newton steps !the condition of the labeling matrix is tested integer bis, i1, i2, k2, count, maxct, k, n, n1, level real stol, kappa, cdmax, newtl, ferr parameter(n = 2, n1 = n+1) !dimension of the problem parameter(ferr = 1.0e-6) !tolerance, used for stopping parameter(stol = 1.e-4) !$\approx$ sqrt(machine tolerance) parameter(cdmax = 1./stol) !maximum for condition estimate parameter(bis = 18) !maximal number of bisections parameter(kappa = 0.5) !contr. factor for Newton steps parameter(maxct = 400) !maximal number of steps real d(0:n1) !level of vertices (stepsize) real v(n, 0:n1) !vertices of simplex integer l(0:n1), r(0:n1) !permutations for vertices integer a(0:n1) !axis from v(.,i) to v(.,i+1) real z(n) !center of virtual simplex real x(0:n), w(0:n) !points , x(0) = 2.**(-level) real x0(n) !starting point real y(n) !current value y= f(x) real c(0:n1) !column real q(0:n1, 0:n1) !orthogonal matrix real b(0:n1, 0:n) !upper triangular matrix logical newl !new level traversed ? logical succ !success of Newton iterations ? i1 = 0 !first higher level index i2 = 0 !last higher level index k2 = 0 !index of vertex being pivoted count = 0 !counts number of function evaluations newtl = 0.5 !last level for Newton steps open(1, file='plhom.dat') !output file call load(v,x0,x,y,z,a,d,r,l,b,q,cdmax,count,n,n1) 1 continue !start of PL loop call index(k2,stol,q,n,n1) !find new pivoting index k2 call pivot(k2,newl,d,v,l,r,a,z,x,i1,i2,n,n1) level = nint(-alog(abs(d(i1)))/alog(2.)) write(1,'(i6,4x,''level='',i3)') count, level write(*,'(i6,4x,''level='',i3)') count, level if (newl .and. (newtl .gt. x(0))) then newtl = x(0) !tentative Newton steps call newton(x,x0,y,w,c,cdmax,ferr,kappa,count, * v,k2,q,b,n,n1,succ) if (succ) then write(1,'(6x,a)') 'Newton iterations succeeded' write(*,'(6x,a)') 'Newton iterations succeeded' goto 35 else write(1,'(6x,a)') 'Newton iterations did not succeed' write(*,'(6x,a)') 'Newton iterations did not succeed' endif do 91 k = 1, n x(k) = v(k, k2) 91 continue endif if (level .gt. bis) then write(1,'(6x,a)') 'maximal bisection level exceeded' write(*,'(6x,a)') 'maximal bisection level exceeded' goto 35 endif call label(x, x0, y, n) count = count + 1 if (count .gt. maxct) then write(1,'(6x,a)') 'maximal number of PL steps exceeded' write(*,'(6x,a)') 'maximal number of PL steps exceeded' goto 35 endif call update(b,q,y,w,cdmax,k2,n,n1) goto 1 !end of PL loop 35 continue !best solution found write(1,'(6x,a,i6/)') 'number of label evaluations:',count write(*,'(6x,a,i6/)') 'number of label evaluations:',count write(1,'(6x,a)') 'approximate solution found:' write(*,'(6x,a)') 'approximate solution found:' write(1,'(6x,a,i2,a,e16.8)') ('x(', k, ')=', x(k), k=1,n) write(*,'(6x,a,i2,a,e16.8)') ('x(', k, ')=', x(k), k=1,n) end subroutine stpnt(x0, n) !user defined !output: x = starting point for homotopy method real x0(n) integer n x0(1) = -1.2 x0(2) = 1.0 return end subroutine stsim(v, n, n1) !user defined !output: v = starting simplex real v(n, 0:n1) integer n, n1, k, m do 91 k = 1, n v(k, 1) = 1.0 91 continue do 92 m = 2, n1 do 93 k = 1, n v(k, m) = v(k, m - 1) 93 continue v(m - 1, m) = 0.0 92 continue return end subroutine label(x, x0, y, n) !user defined !input: x output: y = label of x real x(0:n), x0(n), y(n), x12, x13 integer n, k, level level = nint(-alog(x(0))/alog(2.)) if (level .gt. 0) then !label = f (interesting level) x12 = x(1) * x(1) x13 = x12 * x(1) y(1) = -600.0 * (x(2) - x13) * x12 - 2.0 * (1.0 - x(1)) y(2) = 200.0 * (x(2) - x13) else !label on the trivial level do 91 k = 1, n y(k) = x(k) - x0(k) 91 continue endif return end subroutine givens(b, q, c1, c2, l1, l2, l3, n, n1) !input: b, q, c1, c2, l1, l2, l3 !output: b, q, c1, c2 !one Givens rotation is performed --- !on rows l1 and l2 of b and q !the rotation maps c1, c2 onto sqrt(c1**2+c2**2), 0 real b(0:n1, 0:n),q(0:n1,0:n1),sn,s1,s2,c1,c2,sv1,sv2 integer l1, l2, l3, n, n1, k if (abs(c1)+abs(c2) .eq. 0.) return if (abs(c2) .ge. abs(c1)) then sn = sqrt(1. + (c1/c2)**2) * abs(c2) else sn = sqrt(1. + (c2/c1)**2) * abs(c1) endif s1 = c1/sn s2 = c2/sn do 91 k = 0, n1 sv1 = q(l1, k) sv2 = q(l2, k) q(l1, k) = s1 * sv1 + s2 * sv2 q(l2, k) = -s2 * sv1 + s1 * sv2 91 continue do 92 k = l3, n sv1 = b(l1, k) sv2 = b(l2, k) b(l1, k) = s1 * sv1 + s2 * sv2 b(l2, k) = -s2 * sv1 + s1 * sv2 92 continue c1 = sn c2 = 0.0 return end subroutine testcd(b, cdmax, n, n1) !test of condition --- !a very coarse estimate real b(0:n1, 0:n), cdmax integer n, n1, i, k do 91 i = 1, n do 92 k = 0, i - 1 if (abs(b(k,i)) .gt. cdmax*abs(b(i, i))) then write(1,'(6x,a)') 'bad cond. estimate' write(*,'(6x,a)') 'bad cond. estimate' stop endif 92 continue 91 continue return end subroutine load(v,x0,x,y,z,a,d,r,l,b,q,cdmax,count,n,n1) real v(n,0:n1), x0(n), y(n), z(n), d(0:n1), q(0:n1,0:n1), * b(0:n1,0:n), x(0:n), cdmax integer a(0:n1), l(0:n1), r(0:n1), n, n1, k, m, count call stsim(v, n, n1) call stpnt(x0, n) do 81 k = 1, n y(k) = 0.0 v(k,0) = 0.5 * (v(k,1) + v(k,n1)) !first new vertex do 82 m = 1, n1 y(k) = y(k) + v(k,m) 82 continue y(k) = y(k) / real(n1) !barycenter of starting simplex 81 continue do 83 k = 1, n !shifting barycenter into x0 do 91 m = 0, n1 v(k, m) = v(k, m) - y(k) + x0(k) 91 continue 83 continue do 92 k = 1, n z(k) = 0.5 !load virtual simplex a(k) = k 92 continue do 93 k = 1, n1 d(k) = -1.0 93 continue d(0) = 0.5 do 94 k = 0, n1 r(k) = k + 1 l(k) = k - 1 94 continue l(0) = n1 r(n1) = 0 do 95 m = 0, n1 !loading b and q b(m, n) = 1.0 x(0) = abs(d(m)) do 96 k = 0, n1 q(k, m) = 0.0 96 continue q(m, m) = 1.0 do 97 k = 1, n x(k) = v(k, m) 97 continue call label(x, x0, y, n) count = count + 1 do 98 k = 1, n b(m, k - 1) = y(k) 98 continue 95 continue do 88 m = 0, n do 89 k = m + 1, n1 call givens(b, q, b(m, m), b(k, m), m, k, m+1, n, n1) 89 continue 88 continue call testcd(b, cdmax, n, n1) return end subroutine index(k2,stol,q,n,n1) !find new pivoting index k2 real q(0:n1, 0:n1), s, stol integer k2, n, n1, k if (q(n1, k2).gt. 0) then do 91 k = 0, n1 q(n1, k) = -q(n1, k) 91 continue endif s = 1.e20 k2 = -1 do 92 k = 0, n1 if (q(n1, k) .gt. stol) then if (q(n, k) / q(n1, k) .lt. s) then s = q(n, k) / q(n1, k) k2 = k endif endif 92 continue if (k2 .eq. -1) then write(1,'(6x,a)') 'instability: no index found' write(*,'(6x,a)') 'instability: no index found' stop endif return end subroutine update(b,q,y,w,cdmax,k2,n,n1) !the decomposition q transpose(A) = b is updated --- !for the case that A(k2,1:n) is replaced by y !see section 16.3 real w(0:n), y(n), q(0:n1, 0:n1), b(0:n1, 0:n), cdmax integer k2, n, n1, k, l do 91 k = 1, n w(k - 1) = y(k) 91 continue w(n) = 1.0 do 92 k = 0, n do 82 l = 0, k w(k) = w(k) - b(l, k) * q(l, k2) 82 continue 92 continue do 93 k = n, 0, -1 call givens(b, q, q(k,k2), q(k+1,k2), k, k+1, k, n, n1) 93 continue do 94 k = 0, n1 !correction of q q(0, k) = 0.0 94 continue q(0, k2) = 1.0 do 95 k = 0, n b(0, k) = b(0, k) + w(k) 95 continue do 96 k = 0, n call givens(b, q, b(k,k), b(k+1, k), k, k+1, k+1, n, n1) 96 continue call testcd(b, cdmax, n, n1) return end function even(a) !test if the real number a --- !is near an even number real a, b logical even b = abs(a) / 2.0 if (b - aint(b) .lt. 0.25) then even = .true. else even = .false. endif return end subroutine pivot(k2,newl,d,v,l,r,a,z,x,i1,i2,n,n1) !performs all necessary updates --- !for pivoting vertex k2 of the current simplex --- !with respect to the triangulation $J^3$ !indicates whether a new level is traversed !performs automatic pivots, see section 13.6 --- !in case that the traversed level has height $\le$ 0.25 real d(0:n1), v(n, 0:n1), z(n), x(0:n), s integer l(0:n1), r(0:n1), a(0:n1), k2, i1, i2, n, n1, k1, * k3, k, i0 integer pivcase !cases for pivoting, ordering as in (13.3.7) logical newl newl = .false. 77 continue !entry point for automatic pivot k1 = l(k2) k3 = r(k2) if (k2 .eq. i1) then if (i1 .ne. i2) then pivcase = 2 else pivcase = 5 endif else if (k2 .eq. i2) then if (even((z(a(k1)) + d(k1)) / d(k3))) then pivcase = 7 else pivcase = 8 endif else if (k1 .eq. i2) then if (k3 .ne. i1) then pivcase = 3 else pivcase = 6 endif else if (k3 .eq. i1) then pivcase = 4 else pivcase = 1 endif endif endif endif goto (1,2,3,4,5,6,7,8) pivcase 1 k = a(k2) a(k2) = a(k1) a(k1) = k s = d(k2) d(k2) = d(k1) d(k1) = s do 91 k = 1, n v(k, k2) = v(k, k1) + v(k, k3) - v(k, k2) 91 continue goto 66 2 z(a(k2)) = z(a(k2)) + 2.0 * d(k2) d(k2) = -d(k2) do 92 k = 1, n v(k, k2) = 2.0 * v(k, k3) - v(k, k2) 92 continue goto 66 3 i2 = k2 d(k1) = d(k2) * 0.5 d(k2) = d(k1) a(k1) = a(k2) do 93 k = 1, n v(k, k2) = v(k, k1) + 0.5 * (v(k, k3) - v(k, k2)) 93 continue if ((k3 .eq. l(i1)) .and. (abs(d(k2)) .le. 0.25)) then r(k1) = k3 !automatic pivot l(k3) = k1 r(k3) = k2 l(k2) = k3 r(k2) = i1 l(i1) = k2 d(k2) = d(k3) d(k3) = d(i1) i2 = k3 goto 77 endif goto 66 4 a(i2) = a(k1) d(i2) = -d(k1) * 0.5 d(k2) = d(i2) do 97 k = 1, n v(k, k2) = v(k, i2) + 0.5 * (v(k, k1) - v(k, k2)) 97 continue i3 = r(i2) r(k2) = i3 l(i3) = k2 r(i2) = k2 l(k2) = i2 r(k1) = k3 l(k3) = k1 i2 = k2 if ((r(k2) .eq. k1) .and.(abs(d(k2)) .le. 0.25)) then i2 = l(k2) !automatic pivot r(i2) = k1 l(k1) = i2 r(k1) = k2 l(k2) = k1 r(k2) = k3 l(k3) = k2 i2 = k1 d(k2) = d(k1) d(k1) = d(i1) goto 77 endif goto 66 5 i1 = k3 i2 = k1 d(k2) = d(k2) * 4.0 do 87 k = 1, n v(k, k2) = v(k, k1) 87 continue i0 = l(i1) do 94 k = 0, n1 if ((k .ne. i2) .and. (k .ne. i0)) * z(a(k)) = z(a(k)) - 0.5 * d(k) 94 continue if (abs(d(k2)) .le. 0.5) then !automatic pivot i2 = l(k1) r(i2) = k2 l(k2) = i2 r(k2) = k1 l(k1) = k2 r(k1) = k3 l(k3) = k1 i2 = k2 d(k1) = d(k2) d(k2) = d(k3) goto 77 endif goto 66 6 i1 = k2 i2 = k2 d(k2) = d(k2) * 0.25 newl = .true. do 95 k = 1, n v(k, k2) = 0.5 * (v(k, k1) + v(k, k3)) 95 continue i0 = l(i1) do 96 k = 0, n1 if ((k .ne. i2) .and. (k .ne. i0)) * z(a(k)) = z(a(k)) + 0.5 * d(k) 96 continue goto 66 7 a(k2) = a(k1) d(k2) = d(k1) * 2.0 i2 = k1 do 98 k = 1, n v(k, k2) = v(k, k3) + 2.0 * (v(k, k1) - v(k, k2)) 98 continue goto 66 8 r(k1) = k3 l(k3) = k1 k3 = l(i1) do 99 k = 1, n v(k, k2) = v(k, k3) + 2.0 * (v(k, k1) - v(k, k2)) 99 continue r(k3) = k2 l(k2) = k3 r(k2) = i1 l(i1) = k2 d(k3) = -d(k1) * 2.0 d(k2) = d(k3) a(k3) = a(k1) i2 = k1 goto 66 66 continue !end of pivoting cases do 89 k = 1, n x(k) = v(k, k2) 89 continue x(0) = abs(d(k2)) return end subroutine newton(x,x0,y,w,c,cdmax,ferr,kappa,count, * v,k2,q,b,n,n1,succ) !tentative Newton steps w.r.t. barycentric co-ordinates !see section 13.5 real v(n,0:n1),x(0:n),y(n),q(0:n1,0:n1),b(0:n1,0:n),c(0:n1), * x0(n),w(0:n), cdmax, s, y1, y2, ferr, kappa integer count, k2, n, n1, l, k logical succ succ = .false. s = q(n, k2) / q(n1, k2) !first Newton step do 91 l = 0, n1 c(l) = (q(n, l) - s * q(n1, l)) / b(n, n) 91 continue do 92 k = 1, n x(k) = 0.0 do 93 l = 0, n1 if (l .ne. k2) x(k) = x(k) + c(l) * v(k, l) 93 continue 92 continue call label(x, x0, y, n) count = count + 1 y2 = 0.0 do 94 k = 1, n y2 = y2 + abs(y(k)) 94 continue write(1,'(6x,a,e16.8)') 'norm(f)=', y2 write(*,'(6x,a,e16.8)') 'norm(f)=', y2 77 continue !begin loop of successive Newton steps call update(b,q,y,w,cdmax,k2,n,n1) y1 = y2 s = (1.0 - q(n, k2) / b(n, n)) / q(n1, k2) do 96 l = 0, n1 c(l) = (q(n, l) / b(n, n) + s * q(n1, l)) 96 continue do 97 k = 1, n do 98 l = 0, n1 if (l .ne. k2) x(k) = x(k) + c(l) * v(k, l) 98 continue 97 continue call label(x, x0, y, n) count = count + 1 y2 = 0.0 do 99 k = 1, n y2 = y2 + abs(y(k)) 99 continue write(1,'(6x,a,e16.8)') 'norm(f)=', y2 write(*,'(6x,a,e16.8)') 'norm(f)=', y2 if (y2 .lt. ferr) then succ = .true. return endif if (y2 .gt. kappa * y1) then return else goto 77 endif end C Program 3 program contup !continuation method !follows a curve H(u) = 0 !one Euler predictor, one Newton-corrector !Broyden update after each step, see chapter 7 !stops at a point x such that x(n1) = 0 !arrays: parameter(n = 10, n1 = n+1) !dimension of the problem parameter(pi = 3.1415926535898) dimension b(n1,n) !transpose of Jacobian dimension q(n1,n1) !orth. matrix for QR dec. of b dimension x(n1), u(n1), v(n1) !current points on the curve dimension t(n1) !tangent vector dimension y(n),w(n),p(n),pv(n),r(n) !values of the map H logical test, succ, newton !parameters: ctmax = .8 !maximal contr. rate in corrector step dmax = .2 !maximal norm for H dmin = .001 !minimal norm for H pert = .00001 !perturbation of H hmax = 1.28 !maximal stepsize hmin = .000001 !minimal stepsize hmn = .00001 !minimal Newton step size h = .32 !initial stepsize cdmax = 1000. !maximum for condition estimate angmax = pi/3. !maximal angle maxstp = 9000 !maximal number of evaluations of H acfac = 2. !acceleration factor for steplength control !main program open(1, file='contup.dat') !output file call stpnt(x, n1) !user defined starting point, H(x) = 0 newton = .false. mapct = 0 !counts the calls of the map H call jac(b, x, y, h, n, n1) !b = H$'$(x)$^*$ mapct = mapct + 1 + n1 call decomp(b, q, cond, n, n1) !b, q := orthog. decomp. of b if (cond .gt. cdmax) then write(1,*) ' bad cond. estimate in init. point = ', cond write(*,*) ' bad cond. estimate in init. point = ', cond stop endif do 90 k = 1, n1 !tangent saved t(k) = q(n1, k) 90 continue call setor(or, t, n1) !set orientation 12 continue !begin PC loop if (abs(h).lt.hmin) then write(1,*) ' failure at minimal stepsize' write(*,*) ' failure at minimal stepsize' stop endif if (mapct .gt. maxstp) then write(*,*) ' maximal number of function eval. exceeded' write(1,*) ' maximal number of function eval. exceeded' stop endif do 83 k = 1, n1 !tangent saved t(k) = q(n1, k) 83 continue do 92 k = 1, n1 u(k) = x(k) + h * or * t(k) !predictor step 92 continue call map(u, w, n, n1) mapct = mapct + 1 call upd(q,b,x,u,y,w,t,h,angmax,test,n,n1) !predictor update if (test .eq. .false.) goto 21 !angle test is neg. call newt(q,b,u,v,w,p,pv,r,pert,dmax,dmin, * ctmax,cdmax,test,n,n1) !Newton corrector and update mapct = mapct + 1 if (test.eq..false.) goto 21 !residual or contr. test is neg. goto 22 21 h = h / acfac !PC not accepted goto 12 22 continue !PC step accepted succ = .false. if (v(n1).ge.1.) newton = .true. !switch to Newton steplength if (newton) then h = - (v(n1) - 1.) / q(n1, n1) if (abs(h).lt.hmn) succ = .true. !solution point found else h = abs(h) * acfac !steplength adaptation if (h.gt.hmax) h = hmax endif do 94 k = 1, n1 x(k) = v(k) !new point on curve 94 continue do 95 k = 1, n y(k) = r(k) !y = H(x) 95 continue if (succ) then !stopping the curve tracing write(1,*) ' success with', mapct,' calls of ''map''' write(*,*) ' success with', mapct,' calls of ''map''' write(1,*) write(*,*) write(1,*) ' solution vector:' write(*,*) ' solution vector:' write(1,*) ' ===============' write(*,*) ' ===============' do 96 k = 1, n write(1,*) ' x(', k, ') = ', x(k) write(*,*) ' x(', k, ') = ', x(k) 96 continue stop endif goto 12 end subroutine map(x, y, n, n1) !user defined !input: x output: y = H(x) !H(x) = 0 defines the curve to be traced dimension x(n1), y(n) s = 0. do 91 i = 1, n s = s + x(i) 91 continue do 92 i = 1, n y(i) = x(i) - x(n1) * exp(cos(i * s)) 92 continue return end subroutine jac(b, x, y, h, n, n1) !input: x output: b !evaluates the transpose b of the Jacobian at x !by using forward differences dimension b(n1,n), x(n1), y(n) do 91 i = 1, n1 x(i) = x(i) + h call map(x, y, n, n1) x(i) = x(i) - h do 92 k = 1, n b(i,k) = y(k) 92 continue 91 continue call map(x, y, n, n1) do 93 i = 1, n1 do 94 k = 1, n b(i,k) = (b(i,k) - y(k)) / h 94 continue 93 continue return end subroutine stpnt(x, n1) !user defined !output: x = starting point on curve dimension x(n1) do 91 k = 1, n1 x(k) = 0. 91 continue return end subroutine setor(or, t, n1) !user defined !input: t output: or(t) !decides in which direction the curve will be traversed dimension t(n1) if (t(n1).gt.0.) then or = 1.0 else or = -1.0 endif return end subroutine givens(b, q, c1, c2, l1, l2, l3, n, n1) !input: b, q, c1, c2, l1, l2, l3 !output: b, q, c1, c2 !one Givens rotation is performed --- !on rows l1 and l2 of b and q !the rotation maps c1, c2 onto sqrt(c1**2+c2**2), 0 dimension b(n1, n), q(n1, n1) if (abs(c1)+abs(c2) .eq. 0.) return if (abs(c2) .ge. abs(c1)) then sn = sqrt(1. + (c1/c2)**2) * abs(c2) else sn = sqrt(1. + (c2/c1)**2) * abs(c1) endif s1 = c1/sn s2 = c2/sn do 91 k = 1, n1 sv1 = q(l1, k) sv2 = q(l2, k) q(l1, k) = s1 * sv1 + s2 * sv2 q(l2, k) = -s2 * sv1 + s1 * sv2 91 continue do 92 k = l3, n sv1 = b(l1, k) sv2 = b(l2, k) b(l1, k) = s1 * sv1 + s2 * sv2 b(l2, k) = -s2 * sv1 + s1 * sv2 92 continue c1 = sn c2 = 0. return end subroutine decomp(b, q, cond, n, n1) !input: b output: b, q, cond !a QR decomposition for b is stored in q, b --- !by using Givens rotations on b and q = id --- !until b is upper triangular !a very coarse condition estimate cond is provided dimension b(n1, n), q(n1, n1) do 91 k = 1, n1 !start with q := id do 92 l = 1, n1 q(k, l) = 0. 92 continue q(k, k) = 1.0 91 continue do 93 m = 1, n !successive Givens transformations do 94 k = m+1, n1 call givens(b, q, b(m, m), b(k, m), m, k, m+1, n, n1) 94 continue 93 continue cond = 0. !very coarse condition estimate do 95 i = 2, n do 96 k = 1, i - 1 cond = max(cond, abs(b(k,i)/b(i,i))) 96 continue 95 continue return end subroutine newt(q,b,u,v,w,p,pv,r,pert,dmax,dmin, * ctmax,cdmax,test,n,n1) !input q, b, u, w = H(u) !output v, test, r = H(v) !w is changed !one Newton step v := u - A$^+$ w is performed !where A $\approx$ H$'$ !q, b = QR decomposition of A$^*$ !q, b are updated !perturbations are used for stabilization !residual and contraction tests are performed dimension q(n1,n1),b(n1,n),u(n1),v(n1),w(n),pv(n),p(n),r(n) logical test test = .true. do 81 k = 1, n !perturbation if (abs(w(k)) .gt. pert) then pv(k) = 0. else if (w(k) .gt. 0.) then pv(k) = w(k) - pert else pv(k) = w(k) + pert endif w(k) = w(k) - pv(k) 81 continue d1 = ynorm(w, n) if (d1 .gt. dmax) then test = .false. return endif do 91 k = 1, n do 92 l = 1, k-1 w(k) = w(k) - b(l, k) * w(l) 92 continue w(k) = w(k) / b(k, k) 91 continue d2 = ynorm(w, n) do 93 k = 1, n1 s = 0. do 94 l = 1, n s = s + q(l, k) * w(l) 94 continue v(k) = u(k) - s 93 continue call map(v, r, n, n1) do 74 k = 1, n p(k) = r(k) - pv(k) 74 continue d3 = ynorm(p, n) contr = d3 / (d1 + dmin) if (contr .gt. ctmax) test = .false. do 95 k = n-1, 1, -1 call givens(b, q, w(k), w(k+1), k, k+1, k, n, n1) 95 continue do 96 k = 1, n b(1,k) = b(1,k) - p(k) / d2 96 continue do 97 k = 1, n-1 call givens(b, q, b(k,k), b(k+1,k), k, k+1, k, n, n1) 97 continue if (b(n,n) .lt. 0.) then test = .false. b(n,n) = - b(n,n) do 82 k = 1, n1 q(n,k) = - q(n,k) q(n1,k) = - q(n1,k) 82 continue endif do 85 i = 2, n !perturbation of upper triangular matrix do 86 k = 1, i - 1 if (abs(b(k,i)) .gt. cdmax * abs(b(i,i))) then if (b(i,i) .gt. 0.) then b(i,i) = abs(b(k,i)) / cdmax else b(i,i) = - abs(b(k,i)) / cdmax endif endif 86 continue 85 continue do 87 k = 1, n-1 b(k+1,k) = 0. 87 continue return end subroutine upd(q,b,x,u,y,w,t,h,angmax,test,n,n1) !input q, b, x, u = predictor, y = H(x), w = H(u) !q, b = QR decomposition of transpose(H$'$) !q, b are updated !perturbations are used for stabilization !an angle test is performed dimension q(n1,n1),b(n1,n),x(n1),u(n1),t(n1),y(n),w(n) logical test test = .true. pi = 3.14159265358979323846 do 91 k = 1, n b(n1,k) = (w(k) - y(k)) / h 91 continue do 92 k = 1, n !update call givens(b, q, b(k,k), b(n1,k), k, n1, k, n, n1) 92 continue ang = 0. do 93 k = 1, n1 !angle ang = ang + t(k) * q(n1, k) 93 continue if (ang .gt. 1.0) ang = 1. if (ang .lt. -1.0) ang = -1. ang = acos(ang) if (ang .gt. angmax) test = .false. return end function ynorm(y,n) dimension y(n) s = 0. do 13 k = 1, n s = s + y(k)**2 13 continue ynorm = sqrt(s) return end C Program 4 program bif !continuation method, follows a curve H(u) = 0 !Euler predictor, Newton-correctors !stepsize control by asymptotic estimates !Jacobian is evaluated at each point !an interactive driver monitors: !orientation, perturbation, arclength countdown !in order to trace bifurcating curves !a protocol is written on bif.dat !arrays: parameter(n = 10, n1 = n+1) !dimension of the problem parameter(pi = 3.1415926535898) dimension b(n1,n) !transpose of Jacobian dimension q(n1,n1) !orth. matrix for QR dec. of b dimension x(n1), u(n1) !current points on the curve dimension t(n1) !tangent vector dimension y(n), pv(n) !stores values y := H(x) logical pert, corr, succ !parameters: tol = .0001 !tolerance for corrector iteration ctmax = 0.3 !maximal contr. rate in corrector step dmax = .05 !maximal distance to curve amax = pi/180.*30. !maximal angle hmax = 1. !maximal stepsize hmin = .0001 !minimal stepsize cdmax = 1000. !maximum for condition estimate fmax = 2. !maximal factor for acceleration h = .03 !initial stepsize or = -1. !initial orientation pert = .false. !initial perturbation arc = 0. !initial arclength countdown !main program open(1, file='bif.dat') !output file call stpnt(x, n1) !user defined starting point, H(x) = 0 call setper(pv, n) !set perturbation vectors save = dmax/2./xnorm(pv, n) do 23 k = 1, n !adapt perturbation to dmax pv(k) = pv(k)*save 23 continue 75 call driver(arc,or,pert,h) corr = .true. 12 continue !begin PC loop if (abs(h).lt.hmin) then write(*,*) 'minimal stepsize' write(1,*) 'minimal stepsize' goto 75 endif if (arc.le.0.) then write(*,*) 'enter new arclength for countdown' write(1,*) 'enter new arclength for countdown' goto 75 endif if (corr) then !initial corrector necessary do 82 k = 1, n1 u(k) = x(k) 82 continue call crloop(u,y,pv,b,q,fac, * tol,fmax,ctmax,dmax,cdmax,pert,succ,n,n1) if (succ) then !corrector loop successful? corr = .false. do 18 k = 1, n1 x(k) = u(k) !new point t(k) = q(n1,k) !new tangent 18 continue write(*,*) '||x||=',xnorm(x,n-1),' lambda=',x(n), * ' period=', x(n1), ' h=',h write(1,*) '||x||=',xnorm(x,n-1),' lambda=',x(n), * ' period=', x(n1), ' h=',h else write(*,*) 'initial corrector loop not successful' write(1,*) 'initial corrector loop not successful' goto 75 endif endif do 97 k = 1, n1 u(k) = x(k) + h * or * t(k) !predictor step 97 continue call crloop(u,y,pv,b,q,fac, * tol,fmax,ctmax,dmax,cdmax,pert,succ,n,n1) if (.not.succ) goto 21 angle = 0. do 24 k = 1, n1 angle = angle + t(k)*q(n1,k) 24 continue sangle = sign(1.,angle) angle = sangle*angle if (angle.gt.1.) angle = 1. angle = acos(angle) if ((pert).and.(sangle.lt.0.)) goto 21 if (angle.gt.amax) goto 21 !angle test fac = max(fac, angle/amax*fmax) goto 22 21 h = h / fmax !PC not accepted goto 12 22 continue !PC step accepted arc = arc - abs(h) !arclength countdown if (fac.gt.fmax) fac = fmax h = min(abs(h/fac), hmax) !steplength adaptation if (h.gt.hmax) h = hmax do 94 k = 1, n1 x(k) = u(k) !new point on curve t(k) = q(n1, k) !new tangent 94 continue write(*,*) '||x||=',xnorm(x,n-1),' lambda=',x(n), * ' period=', x(n1), ' h=',h write(1,*) '||x||=',xnorm(x,n-1),' lambda=',x(n), * ' period=', x(n1), ' h=',h if ((.not.pert).and.(sangle.lt.0.)) then write(*,*) 'orientation changed' write(1,*) 'orientation changed' goto 75 endif goto 12 end subroutine driver(arc,or,pert,h) !interactive driver logical pert realn = 0. write(1,*) ' ' write(1,*) 'interactive driver' 77 continue write(*,*) ' 1) stop 2) go 3) arc=',arc,' 4) or=',or, * ' 5) pert=',pert,' 6) h=',h write(*,*) 'enter integer (option) and realnumber (value)' read(*,*) intgr, realn write(1,*) 'integer=',intgr, ' real number=', realn if (intgr .eq. 1) stop if (intgr .eq. 2) goto 78 if (intgr .eq. 3) arc = realn if (intgr .eq. 4) or = -or if (intgr .eq. 5) pert = .not.pert if (intgr .eq. 6) h = realn goto 77 78 write(1,*) 'arc=', arc, ' or=', or, ' pert=', pert, ' h=', h write(1,*) ' ' end subroutine map(x, y, n, n1) !user defined !input: x output: y = H(x) !H(x) = 0 defines the curve to be traced dimension x(n1), y(n), w(200), d(200), xm(200), * p(200), q(200) nm3 = n-3 nm2 = n-2 nm1 = n-1 lim = 4*nm2 + n h = 1./ float(nm2) xind = x(n1)/ 2./ h ind = xind if ((ind.lt.1).or.(ind+n.gt.lim)) goto 2 t = (xind - ind)*h r1 = -x(n)*h*5./12. r2 = -x(n)*h*8./12. r3 = x(n)*h/12. q1 = -x(n)*h/3. q2 = -x(n)*h*4./3. q3 = q1 do 1 k = 1, lim if (k.le.nm1) then w(k) = x(k) elseif (mod(k,2).eq.0) then w(k) = w(k-1) + r1* xf(w(k-nm1)) + r2* xf(w(k-nm2)) * + r3* xf(w(k-nm3)) else w(k) = w(k-2) + q1* xf(w(k-n)) + q2* xf(w(k-nm1)) * + q3* xf(w(k-nm2)) endif 1 continue do 3 k = 2, lim - 1 d(k) = 3./h**2 * ( w(k+1)-2.*w(k)+w(k-1) ) 3 continue p(2) = sqrt(2.) d(2) = d(2)/ p(2) do 4 k = 3, lim - 1 q(k-1) = .5/ p(k-1) p(k) = sqrt(2. - q(k-1)**2) d(k) = (d(k) - q(k-1)*d(k-1))/ p(k) 4 continue xm(lim) = 0. xm(lim-1) = d(lim-1)/ p(lim-1) do 5 k = lim - 2, 2, -1 xm(k) = (d(k) - q(k)*xm(k+1))/ p(k) 5 continue xm(1) = 0. do 6 k = ind + 1, ind + nm1 a7 = w(k) c7 = xm(k)/2. b7 = (w(k+1)-w(k))/h - h/6.*(2.*xm(k)+xm(k+1)) d7 = (xm(k+1)-xm(k))/ (6.*h) y(k-ind) = x(k-ind) + (((d7*t)+c7)*t+b7)*t+a7 6 continue y(n) = x(1) return 2 write(*,*) 'failure in map' write(1,*) 'failure in map' stop end function xf(t) !auxiliary function for above map xf = t*(1. + t**2) / (1. + t**4) end subroutine jacob(b, x, n, n1) !user defined !input: x output: b !evaluates the transpose b of the Jacobian at x dimension b(n1,n), x(n1), y(30), w(30) h1 = 1024. h = 1./h1 call map(x, y, n, n1) do 1 k = 1, n1 x(k) = x(k) + h call map(x, w, n, n1) x(k) = x(k) - h do 2 l = 1, n b(k,l) = h1*(w(l)-y(l)) 2 continue 1 continue end subroutine stpnt(x, n1) !user defined !output: x = starting point on curve parameter(pi = 3.1415926535898) dimension x(n1) h = 1./float(n1-3) do 1 k=1, n1-2 tk = (k-1)*h*pi/2. x(k) = .1* sin(tk) 1 continue x(n1-1) = pi/2. x(n1) = 4. end subroutine setper(pv, n) !user defined !defines the perturbation vector pv dimension pv(n) do 1 k = 1, n-1 pv(k) = 1. 1 continue pv(n) = 0. end subroutine func(x, y, pv, pert, n, n1) !perturbed function evaluation dimension x(n1), y(n), pv(n) logical pert call map(x, y, n, n1) if (pert) then do 31 k = 1, n y(k) = y(k) - pv(k) 31 continue endif end subroutine givens(b, q, c1, c2, l1, l2, l3, n, n1) !input: b, q, c1, c2, l1, l2, l3 !output: b, q, c1, c2 !one Givens rotation is performed --- !on rows l1 and l2 of b and q !the rotation maps c1, c2 onto sqrt(c1**2+c2**2), 0 dimension b(n1, n), q(n1, n1) if (abs(c1)+abs(c2) .eq. 0.) return if (abs(c2) .ge. abs(c1)) then sn = sqrt(1. + (c1/c2)**2) * abs(c2) else sn = sqrt(1. + (c2/c1)**2) * abs(c1) endif s1 = c1/sn s2 = c2/sn do 91 k = 1, n1 sv1 = q(l1, k) sv2 = q(l2, k) q(l1, k) = s1 * sv1 + s2 * sv2 q(l2, k) = -s2 * sv1 + s1 * sv2 91 continue do 92 k = l3, n sv1 = b(l1, k) sv2 = b(l2, k) b(l1, k) = s1 * sv1 + s2 * sv2 b(l2, k) = -s2 * sv1 + s1 * sv2 92 continue c1 = sn c2 = 0.0 end subroutine decomp(b, q, cond, n, n1) !input: b, output: b, q, cond !a QR decomposition for b is stored in q, b --- !by using Givens rotations on b and q = id --- !until b is upper triangular !a very coarse condition estimate cond is provided dimension b(n1, n), q(n1, n1) do 91 k = 1, n1 !start with q := id do 92 l = 1, n1 q(k, l) = 0.0 92 continue q(k, k) = 1.0 91 continue do 93 m = 1, n !successive Givens transformations do 94 k = m+1, n1 call givens(b, q, b(m, m), b(k, m), m, k, m+1, n, n1) 94 continue 93 continue cond = 0. !very coarse condition estimate do 95 i = 2, n do 96 k = 1, i - 1 cond = max(cond, abs(b(k,i)/b(i,i))) 96 continue 95 continue end subroutine newton(q, b, u, y, n, n1) !input q, b, u, y = H(u), n, n1 !output u !y is changed !a Newton step $u := u - A^+ (H(u)-pv)$ is performed --- !where A approximates the current Jacobian $H'$ !q, b = QR decomposition of $A^*$ dimension q(n1, n1), b(n1, n), u(n1), y(n) do 91 k = 1, n do 92 l = 1, k-1 y(k) = y(k) - b(l, k) * y(l) 92 continue y(k) = y(k) / b(k, k) 91 continue do 93 k = 1, n1 s = 0.0 do 94 l = 1, n s = s + q(l, k) * y(l) 94 continue u(k) = u(k) - s 93 continue end subroutine crloop(x,y,pv,b,q,fac, * tol,fmax,ctmax,dmax,cdmax,pert,succ,n,n1) !corrector loop !input x,y !output x,y,b,q,fac,succ dimension x(n1),y(n),pv(n),b(n1,n),q(n1,n1) logical succ, pert succ = .false. !success of corrector loop fac = 1./fmax call func(x, y, pv, pert, n, n1) 35 continue !begin loop call jacob(b, x, n, n1) !b := transpose of Jacobian at x call decomp(b, q, cond, n, n1) !decompose b if (cond .gt. cdmax) return !bad conditioning dist1 = xnorm(y,n) fac = max(fac, sqrt(dist1/dmax)*fmax) if (dist1.lt.tol) goto 34 !corrector successful if (dist1.gt.dmax) return call newton(q, b, x, y, n, n1) call func(x, y, pv, pert, n, n1) dist2 = xnorm(y,n) contr = dist2 / (dist1 + tol) !contraction rate fac = max(fac, sqrt(contr/ctmax)*fmax) if (contr.gt.ctmax) return dist1 = dist2 goto 35 !end loop 34 succ = .true. !corrector successful end function xnorm(y, n) !calculates euclidean norm of y dimension y(n) x = 0. do 1 k = 1, n x = x + y(k)**2 1 continue xnorm = sqrt(x) end C Program 5 program surapp!Piecewise linear approximation of an !implicitly defined surface, see section 15.4 parameter (n=1, k=2, nplusk=n+k)!dimension of the problem parameter (k1=k+1, n3=n+3) parameter (lisdim = 1000)!dimension of the simplex list double precision eps parameter (eps = 1d-10)!machine tolerance double precision delta, origin(1:nplusk) !mesh size and origin of the triangulation double precision lbound(1:nplusk), ubound(1:nplusk) !lower and upper bounds of the problem double precision simx(0:nplusk,1:nplusk) !vertices of current simplex double precision simf(0:nplusk,1:n) !function values of vertices double precision error!error of the approximation double precision u(1:nplusk), v(1:nplusk), x(1:nplusk), * fx(1:n)!auxiliary arrays integer slist (1:lisdim,1:nplusk)!list of simplices !a simplex is characterized by barycenter * (nplusk+1) integer inds!current member of simplex list integer maxs!last entry of simplex list integer numver (k1:n3)!counts pl pieces integer pi(1:nplusk), z(1:nplusk) !pi and z values of the current simplex integer clface (0:nplusk), i!auxiliary variables logical facets (0:nplusk) !transverse facets of the current simplex !facets(i) = .true. means that facet i is transverse open(1, file='surapp.dat') !output file inds = 1!starting values (initialization) maxs = 1 error = 0.0 do 10 i = k1, n3 numver (i) = 0 10 continue write (*,'(/1x,a,a)') 'pl approximation of an implicitly', * ' defined surface' write (1,'(/1x,a,a)') 'pl approximation of an implicitly', * ' defined surface' call start (delta, origin, lbound, ubound, simx, simf, * slist, nplusk, n, lisdim, pi, z, x, v, fx, u) !compute starting simplex, mesh size, origin, lower and !upper bounds 20 continue!begin of pl loop call appsim (inds, simx, simf, numver, nplusk, n, k1, n3, * slist, lisdim, maxs, lbound, ubound, eps, * error, pi, z, u, v, x, fx, facets, clface ) !process current simplex inds = inds+1 if ((inds .le. maxs) .and. (inds.le.lisdim) ) then !not all simplices are processed call getver (inds, simx, simf, slist, nplusk, lisdim, * n, pi, z, x, fx, u, origin, delta) !simx and simf of next simplex are computed goto 20 end if!end of loop !statistics of program write (*,'(//'' total number of transverse simplices'', * 14x,i8)') maxs write (*,'(5x,''pl pieces containing'',i2,a,15x,i8)') * (i,' vertices',numver(i),i=k1,n3) write (*,'(/'' maximum of all function values'',16x,d12.6)') * error write (1,'(//'' total number of transverse simplices'', * 14x,i8)') maxs write (1,'(5x,''pl pieces containing'',i2,a,15x,i8)') * (i,' vertices',numver(i),i=k1,n3) write (1,'(/'' maximum of all function values'',16x,d12.6)') * error stop end subroutine appsim (inds, simx, simf, numver, nplusk, n, k1, * n3, slist, lisdim, maxs, lbound, ubound, * eps, error, pi, newcen, u, v, x, fx, * facets, clface ) !input: inds, simx, simf, numver, slist, maxs, lbound, ubound !output: numver, slist, maxs !this subprogram computes all cl faces of the current simplex, !all neighbors of the current simplex which share a common !transverse facet are put on the simplex list double precision simx (0:nplusk,1:nplusk), * simf (0:nplusk,1:n), * lbound (1:nplusk), ubound(1:nplusk), * eps, error, * u(1:nplusk), v(1:nplusk), * x(1:nplusk), fx(1:n) integer i, j, numver (k1:n3), lisdim, k1, n3, maxs, inds, * slist (1:lisdim,1:nplusk), pi(1:nplusk) logical facets (0:nplusk) !for an explanation of these variables see the main program integer clface (0:nplusk)!indices of cl face (0..n) and of !vertices to be pivoted (n+1..n+k) integer start!first vertex to pivot in cl face integer numcl!counts cl faces integer newcen (1:nplusk)!barycenter * (nplusk+1) of a neighbor logical bound!function which checks the bounds !search of a transverse edge !works only if nplusk=3 i = 1 110 continue if ((simf(0,1).le.0.0) .eqv. (simf(i,1).le.0.0)) then i = i+1 if (i .le. nplusk) goto 110 end if if (i .gt. nplusk) return do 120 j = 0, nplusk!starting values (initialization) clface(j) = j 120 continue numcl = 1 if (i.ne.1) then clface(1) = i clface(i) = 1 end if do 130 i = 0, nplusk !facets clface(n+1) and clface(n+2) are transverse facets (clface(i)) = i .gt. n 130 continue start = clface (n+1) call output (simx, simf, clface, u, v, x, fx, inds, * nplusk, n, eps, error) !compute zero point of the cl face !and write it on a file or screen 140 continue!begin cl face loop call pivot (clface, simf, nplusk, n) !compute next cl face call output (simx, simf, clface, u, v, x, fx, inds, * nplusk, n, eps, error) !compute zero point of the cl face !and write it on a file or screen numcl = numcl+1 facets (clface(nplusk)) = .true. !facet clface(nplusk) is transverse if (clface(n+2) .ne. start) goto 140 !stop test works correctly if k=2 numver(numcl) = numver(numcl)+1!counts pl pieces do 160 i = 0, nplusk !loop which checks the bounds of transverse facets !a facet is outside the boundary if all vertices of the !facet are outside if (facets(i)) then do 150 j = 0, nplusk if ((j .ne. i) .and. .not. * bound (simx, j, lbound, ubound, nplusk, eps)) * goto 160 150 continue facets (i) = .false. end if 160 continue do 170 i = 0, nplusk !loop which computes all neighbors of the current simplex !if they share a common transverse facet if (facets(i)) then call reflec(slist,inds,i,newcen,nplusk,lisdim,pi) !compute the barycenter of the neighbor and !put it on the simplex list call newsim (newcen, slist, maxs, lisdim, nplusk) end if 170 continue return end logical function bound(simx,ind,lbound,ubound,nplusk,eps) !input: simx, ind, lbound, ubound !output: bound !this function checks the bounds of the vertex simx(ind,.) double precision simx (0:nplusk,1:nplusk), * lbound(1:nplusk), ubound(1:nplusk), * eps integer nplusk, ind !for an explanation of these variables see the main program integer i!auxiliary variables logical r i = 0 210 i = i+1 r = ((lbound(i)-simx(ind,i)) .ge. eps) .or. * ((simx(ind,i)-ubound(i)) .ge. eps) if ((i .lt. nplusk) .and. .not. r) goto 210 bound = r return end subroutine fvalue (x, f, nplusk, n)!input: x !output: f (function value of x) !user defined function e.g. an equation of a torus double precision x(1:nplusk), f(1:n), help integer nplusk, n help = x(2)**2 + x(3)**2 f(1) = (x(1)**2 + help + 0.96d0)**2 - 4.0d0*help return end subroutine getver(inds,simx,simf,slist,nplusk,lisdim,n, * pi,z,x,fx,vertex,origin,delta) !input: inds, slist, origin, delta !output: simx, simf !the subroutine computes the vertices of the current simplex !and the function values belonging to the vertices !see the rules of the Freudenthal triangulation (12.1.10) double precision simx(0:nplusk,1:nplusk), * simf(0:nplusk,1:n), * x(1:nplusk), fx(1:n), delta, * origin(1:nplusk) integer lisdim, nplusk, n, slist (1:lisdim,1:nplusk), * inds, pi(1:nplusk), z(1:nplusk) !for an explanation of these variables see the main program integer i, help double precision vertex (1:nplusk)!auxiliary variables do 410 i = 1, nplusk !permutation pi and integer vector z are calculated !(only the barycenter * (nplusk+1) was stored) z(i) = slist(inds,i) / (nplusk+1) help = mod (slist(inds,i),nplusk+1) if (help .lt. 0) then help = help+nplusk+1 z(i) = z(i)-1 end if pi (nplusk+1-help) = i 410 continue do 420 i = 1, nplusk!starting value for current vertex vertex(i) = z(i) 420 continue call simtox (vertex, simx, origin, delta, nplusk, 0) !calculate coordinates of vertex 0 and put it on simx(0,.) !function value of vertex 0 is computed and stored in simf(0,.) do 430 i = 1, nplusk x(i) = simx(0,i) 430 continue call fvalue (x, fx, nplusk, n) do 440 i = 1, n simf(0,i) = fx(i) 440 continue !all other vertices and function values are calculated do 470 i = 1, nplusk vertex (pi(i)) = vertex (pi(i)) + 1.0d0 !rules of the Freudenthal triangulation, see (12.1.10) call simtox (vertex, simx, origin, delta, nplusk, i) do 450 j = 1, nplusk x(j) = simx(i,j) 450 continue call fvalue (x, fx, nplusk, n) do 460 j = 1, n simf(i,j) = fx(j) 460 continue 470 continue return end subroutine newsim (center, slist, maxs, lisdim, nplusk) !input: center, slist, smax !output: slist, smax !the subroutine puts a new barycenter on the simplex list integer maxs, lisdim, nplusk, slist (1:lisdim,1:nplusk) !for an explanation of these variables see the main program integer i, j!auxiliary variables integer center (1:nplusk) !barycenter * (nplusk+1) of a simplex !loop compares all simplices of the list with current center do 520 i = maxs, 1, -1 do 510 j = 1, nplusk if (slist(i,j) .ne. center(j)) goto 520 510 continue return!the simplex is already a member of the list 520 continue !center belongs to a new simplex !and must be stored in the simplex list maxs = maxs+1 !check the capacity of the simplex list if (maxs .eq. lisdim+1) then write (*,'('' simplex list is too small'')') write (1,'('' simplex list is too small'')') return end if !storing center at the end of the simplex list do 530 i = 1, nplusk slist (maxs,i) = center(i) 530 continue return end subroutine output (simx, simf, clface, u, v, x, fx, inds, * nplusk, n, eps, error) !input: simx, simf, clface, inds, error !output: error, zero point on the screen !output calculates the zero point on the edge with !a bisection method and writes it to a file or screen !subroutine works correctly if nplusk=3 double precision simx(0:nplusk,1:nplusk), * simf(0:nplusk,1:n), eps, error integer inds, nplusk, n !for an explanation of these variables see the main program integer clface (0:nplusk) !for an explanation of clface see the subroutine appsim double precision u(1:nplusk) !first vertex (simx(clface(0),.)) double precision v(1:nplusk)!u+v = second vertex double precision x(1:nplusk)!zero point approximation double precision fx(1:n)!function value of x double precision lambda !barycentric coordinate of the zero point double precision lowerb, upperb !lower and upper bound for the bisection method logical neg!=.true. iff f(u) is negative integer i!auxiliary variable !starting values are calculated do 610 i = 1, nplusk u(i) = simx (clface(0),i) v(i) = simx (clface(1),i) - u(i) 610 continue neg = simf(clface(0),1) .lt. 0.0d0 lowerb = 0.0d0 upperb = 1.0d0 620 continue!begin loop of bisection method lambda = (lowerb + upperb)/2.0d0 do 630 i = 1, nplusk x(i) = lambda*v(i) + u(i) 630 continue call fvalue ( x, fx, nplusk, n) if (neg .eqv. (fx(1) .lt. 0.0d0)) then lowerb = lambda else upperb = lambda end if if (upperb-lowerb .ge. eps) goto 620 !approximation error if (error .lt. dabs(fx(1))) error = dabs (fx(1)) write (*,'(i6,5(3x,d15.8))') inds,(x(i),i=1,nplusk) write (1,'(i6,5(3x,d15.8))') inds,(x(i),i=1,nplusk) return end subroutine pivot (clface, simf, nplusk, n) !input: clface, simf !output: clface !pivot of clface(n+1) !pivot works correctly if nplusk = 3 double precision simf (0:nplusk,1:n) integer nplusk, n !for an explanation of these variables see the main program integer clface (0:nplusk) !for an explanation of clface see the subroutine appsim integer leave!index of vertex which leaves the cl face if ((simf(clface (n+1),1).le.0.0d0) .eqv. * (simf(clface (0),1).le.0.0d0)) then !sign of clface(0) equal to sign of clface(n+1) leave = clface(0) clface(0) = clface (n+1) else !sign of clface(1) equal to sign of clface(n+1) leave = clface(1) clface(1) = clface (n+1) end if clface (n+1) = clface (n+2) clface (n+2) = leave return end subroutine reflec (slist, inds, facet, newcen, nplusk, * lisdim, pi) !input: slist, inds, facet !output: newcen !newcen is obtained by reflecting the vertex facet !of the current simplex !see rules (12.1.10) of the Freudenthal triangulation integer nplusk, slist(1:lisdim,1:nplusk), pi(1:nplusk) !for an explanation of these variables see the main program integer facet!index of vertex which should be reflected !from the current simplex integer newcen(1:nplusk) !barycenter * (nplusk+1) of the neighbor integer i, help!auxiliary variables !computing of starting values do 810 i = 1, nplusk newcen(i) = slist (inds,i) help = mod (slist(inds,i),nplusk+1) if (help .lt. 0) help = help+nplusk+1 pi (nplusk+1-help) = i 810 continue !reflection (see (12.1.11)) if ((facet.gt.0) .and. (facet.lt.nplusk)) then newcen (pi(facet)) = newcen (pi(facet)) - 1 newcen (pi(facet+1)) = newcen (pi(facet+1)) + 1 else if (facet.eq.0) then newcen (pi(1)) = newcen(pi(1)) + 2 do 820 i = 2, nplusk newcen(pi(i)) = newcen(pi(i)) + 1 820 continue else newcen (pi(nplusk)) = newcen (pi(nplusk)) - 2 do 830 i = 1, nplusk-1 newcen (pi(i)) = newcen (pi(i)) - 1 830 continue end if return end subroutine simtox (vertex, simx, origin, delta, * nplusk, ind) !input: vertex, origin, delta !output: simx(ind,.) !transformation of vertex to true co-ordinates double precision delta, origin(1:nplusk), * simx(0:nplusk,1:nplusk) integer nplusk !for an explanation of these variables see the main program double precision vertex(1:nplusk) !integer coordinates of vertex integer ind!index for simx integer i do 910 i = 1, nplusk simx(ind,i) = origin(i) + delta*vertex(i) 910 continue return end subroutine start (delta, origin, lbound, ubound, simx, * simf, slist, nplusk, n, lisdim, pi, z, * x, v, fx, u) !output: delta, origin, lbound, ubound, simx, simf, slist !user defined subroutine !calculating starting values of the algorithm double precision delta, origin(1:nplusk), * lbound(1:nplusk), ubound(1:nplusk), * simx(0:nplusk,1:nplusk), * simf(0:nplusk,1:n), x(1:nplusk), * v(1:nplusk), fx(1:n), u(1:nplusk) integer lisdim, nplusk, n, slist (1:lisdim,1:nplusk), * pi(1:nplusk), z(1:nplusk) !for an explanation of these variables see the main program integer i, j, step, max!auxiliary variables !setting bounds lbound(1) = -0.6d0 lbound(2) = -1.2d0 lbound(3) = -1.2d0 ubound(1) = 0.6d0 ubound(2) = 1.2d0 ubound(3) = 1.2d0 !initial point (should be a approximation of a zero point) x(1) = 0.1875d0 x(2) = 1.125d0 x(3) = 0.0625d0 delta = 0.25d0!mesh size write (*,'(/'' initial point'')') write (*,'(5(1x,d10.5))') (x(i),i=1,nplusk) write (1,'(/'' initial point'')') write (1,'(5(1x,d10.5))') (x(i),i=1,nplusk) !first simplex: z = 0 and pi = id do 1010 i = 1, nplusk slist(1,i) = nplusk+1-i 1010 continue step = 0 max = 1 !construction of of a transverse simplex 1020 continue do 1030 i = 1, nplusk!Freudenthal triangulation origin(i) = x(i) - delta/(nplusk+1)*(nplusk+1-i) 1030 continue call getver (1, simx, simf, slist, nplusk, lisdim, n, * pi, z, v, fx, u, origin, delta) !simx and simf of the starting simplex are calculated !search for a transverse edge !works only if nplusk=3 i = 1 1040 continue if ((simf(0,1).le.0.0) .eqv. (simf(i,1).le.0.0)) then i = i+1 if (i .le. nplusk) goto 1040 end if if (i .gt. nplusk) then!simplex is not transverse step = step+1 if (step .lt. max) then !reduce mesh size and try it again delta = delta*step/(step+1.0) goto 1020 else stop end if end if !recording simplex and function values write (*,'(/'' start simplex of mesh size '',f10.5)') * delta write (1,'(/'' start simplex of mesh size '',f10.5)') * delta do 1050 j = 1, nplusk write (*,'(6(1x,f11.5))') (simx (i,j), i=0,nplusk) write (1,'(6(1x,f11.5))') (simx (i,j), i=0,nplusk) 1050 continue write (*,'(/'' function values'')') write (1,'(/'' function values'')') do 1060 j = 1, n write (*,'(6(1x,d11.5))') (simf (i,j), i=0,nplusk) write (1,'(6(1x,d11.5))') (simf (i,j), i=0,nplusk) 1060 continue write (*,'(/'' simplex numbers and approximate zero '', * ''points:'')') write (1,'(/'' simplex numbers and approximate zero '', * ''points:'')') return end