subroutine evupd(n,i,d,z,delta,rho,dlam,ifail) CVD\$G NOCONCUR c*********************************************************************** c this subroutine computes the updated eigenvalues of a c rank one modification to a symmetric matrix. it is assumed c that the eigenvalues are in the array d, and that c c d(i) .ne. d(j) for i .ne. j c c it is also assumed that the eigenvalues are in increasing c order and that the value of rho is positive. this is c arranged by the calling subroutine sesupd, and is no loss c in generality. c it is also assumed that the values in the array z are c the squares of the components of the updatingvector. c c c the method consists of approximating the rational functions c c rho = sum(z(j)/((d(j)-d(i))/rho - lamda): j = i+1,n) c c phi =sum(z(j)/(d(j)-d(i))/rho - lamda): j = 1,i) c c by simple interpolating rational functions. this avoids c the need for safeguarding by bisection since c the convergence is monotone, and quadratic from any starting c point that is greater than zero but less than the solution c c c input variables... c n the length of all arrays c c i the i - th eigenvalue is computed c c d the original eigenvalues. it is assumed that they are c in order, d(i) .lt. d(j) for i .lt. j. c c z this array of length n contains the squares of c of the components of the updating vector c c delta this array of length n contains (d(j) - lamda(i))/rho c in its j - th component. these values will be used c to update the eigenvectors in sesupd. c c rho this is the scalar in the the symmetric updating c formula. c c dlam this scalar will contain the value of the i - th c updated eigenvalue on return c c ifail this integer variable indicates failure of the c updating process with value 1, and success c with value 0. c c c c************************************************************************** integer i,n,im1,ip1,ip2,niter dimension d(n),z(n), delta(n) real d,z,rho,delta,zero,one,two, 1 del,phi,dphi,psi,dpsi,lambda,oldlam, 2 t,temp,a,b,d1,w,eps,eps1,eta,dlam,dmax zero = 0.0e0 one = 1.0e0 two = 2.0e0 c c eps is machine precision c eta is the relative accuracy requirement on the roots. c c these values should be adjusted for the particular machine. c eps = 100.*epslon(1.0) eta = eps*8. c im1 = i - 1 ip1 = i + 1 ip2 = i + 2 niter = 1 lambda = zero oldlam = zero del = d(i) do 100 j = 1,n delta(j) = (d(j) - del)/rho 100 continue dmax = amax1(abs(d(1)), abs(d(n))) dmax = abs(rho) + dmax c c calculate initial guess c if (i .lt. n) * then del = d(ip1) else del = d(n) + rho endif a = zero do 200 j = 1,im1 a = a + rho*z(j)/(d(j) - del) 200 continue b = zero do 220 j = ip2,n b = b + rho*z(j)/(d(j) - del) 220 continue a = a + (b + one) if (ip1 .gt. n) * then t = z(i)/a else t = a*delta(ip1) b = t + z(i) + z(ip1) if (b .ge. zero) * then t = two*z(i)*delta(ip1)/ * (b + sqrt(abs(b*b - 4.0e0*t*z(i)))) else t = (b - sqrt(b*b - 4.0e0*t*z(i)))/(two*a) endif endif c if (t .le. zero) write(6,*) ' int ges ',t c c test to see that the initial guess is not too close to endpoint c if (ip1 .le. n .and. t .ge. .9*delta(ip1)) * t = .9*delta(ip1) c c update the values of the array delta c 250 continue tsave = abs(t) do 300 j = 1,n delta(j) = delta(j) - t 300 continue lambda = lambda + t dlam = d(i) + rho*lambda c c evaluate psi and the derivative dpsi c dpsi = zero psi = zero do 400 j = 1,i t = z(j)/delta(j) psi = psi + t dpsi = dpsi + t/delta(j) 400 continue c c evaluate phi and the derivative dphi c dphi = zero phi = zero if (i .eq. n) go to 600 do 500 j = ip1,n t = z(j)/delta(j) phi = phi + t dphi = dphi + t/delta(j) 500 continue c c test for convergence c return if the test is satisfied c 600 continue w = one + phi + psi eps1 = eps*dmax*sqrt(dphi + dpsi) if ((abs(w) .le. eps1) .and. (abs(lambda - oldlam) .le. 1 eta*oldlam) ) then return endif c c return with ifail = -1 if convergence has not ocurred c within 45 iterations. c if (niter .lt. 45) go to 650 ifail = -1 return 650 continue niter = niter + 1 oldlam = lambda c c calculate the new step c if (i .ne. n) go to 700 c c otherwise c calculate the step for the special case i = n c t = (w*psi)/dpsi go to 250 700 continue del = delta(ip1) temp = psi/dpsi d1 = one + phi - del*dphi a = (del*(one + phi) + psi*temp)/d1 + temp b = (two*temp*del*w)/d1 t = sqrt(abs(a*a - two*b)) t = b/(a + t) if (t .lt. -eps) t = -tsave/two go to 250 c c last card of evupd c end subroutine rotate(icount,indx,jm1,n,ldq,q,z,rhoge0) CVD\$G NOCONCUR c c rotate z components to 0 c assumes norm(z) .eq. 1 c logical rhoge0 integer icount,indx(*),jm1,n,ldq real q(ldq,*),z(*),gamma,sigma,tau,t1,t2 c if (rhoge0) *then do 200 j = jm1-icount+1,jm1-1 jp1 = j + 1 gamma = z(jp1) sigma = z(j) tau = sqrt(gamma*gamma + sigma*sigma) if (tau .gt. 0.0e0) * then gamma = gamma/tau sigma = sigma/tau z(jp1) = tau z(j) = 0.0 jj = indx(j) jjp1 = indx(jp1) do 100 i = 1,n t1 = q(i,jj) t2 = q(i,jjp1) q(i,jj) = t1*gamma - t2*sigma q(i,jjp1) = t1*sigma + t2*gamma 100 continue endif 200 continue else do 400 j = jm1-1,jm1-icount+1,-1 jp1 = j + 1 gamma = z(j) sigma = z(jp1) tau = sqrt(gamma*gamma + sigma*sigma) if (tau .gt. 0.0e0) * then gamma = gamma/tau sigma = sigma/tau z(j) = tau z(jp1) = 0.0e0 jj = indx(j) jjp1 = indx(jp1) do 300 i = 1,n t1 = q(i,jj) t2 = q(i,jjp1) q(i,jjp1) = t2*gamma - t1*sigma q(i,jj) = t2*sigma + t1*gamma 300 continue endif 400 continue endif c return c c last card of rotate c end function rand(foo) CVD\$G NOCONCUR data init/1365/ init = mod(3125*init,65536) rand = (init - 32768.)/(2.**15) return end real function enorm(n,x) CVD\$G NOCONCUR integer n real x(n) tmax = 0.0 do 100 j = 1,n temp = max(tmax,abs(x(j))) if (temp .gt. tmax) tmax = temp 100 continue if (tmax .eq. 0.0) then enorm = tmax else temp = 0.0 do 200 j = 1,n temp = temp + (x(j)/tmax)**2 200 continue enorm = tmax*sqrt(temp) endif return c end subroutine tql2(nm,n,d,e,z,ierr) CVD\$G NOCONCUR c integer i,j,k,l,m,n,ii,l1,l2,nm,mml,ierr real d(n),e(n),z(nm,n) real c,c2,c3,dl1,el1,f,g,h,p,r,s,s2,tst1,tst2,pythag c c this subroutine is a translation of the algol procedure tql2, c num. math. 11, 293-306(1968) by bowdler, martin, reinsch, and c wilkinson. c handbook for auto. comp., vol.ii-linear algebra, 227-240(1971). c c this subroutine finds the eigenvalues and eigenvectors c of a symmetric tridiagonal matrix by the ql method. c the eigenvectors of a full symmetric matrix can also c be found if tred2 has been used to reduce this c full matrix to tridiagonal form. c c on input c c nm must be set to the row dimension of two-dimensional c array parameters as declared in the calling program c dimension statement. c c n is the order of the matrix. c c d contains the diagonal elements of the input matrix. c c e contains the subdiagonal elements of the input matrix c in its last n-1 positions. e(1) is arbitrary. c c z contains the transformation matrix produced in the c reduction by tred2, if performed. if the eigenvectors c of the tridiagonal matrix are desired, z must contain c the identity matrix. c c on output c c d contains the eigenvalues in ascending order. if an c error exit is made, the eigenvalues are correct but c unordered for indices 1,2,...,ierr-1. c c e has been destroyed. c c z contains orthonormal eigenvectors of the symmetric c tridiagonal (or full) matrix. if an error exit is made, c z contains the eigenvectors associated with the stored c eigenvalues. c c ierr is set to c zero for normal return, c j if the j-th eigenvalue has not been c determined after 30 iterations. c c calls pythag for sqrt(a*a + b*b) . c c questions and comments should be directed to burton s. garbow, c mathematics and computer science div, argonne national laboratory c c this version dated august 1983. c c ------------------------------------------------------------------ c ierr = 0 if (n .eq. 1) go to 1001 c do 100 i = 2, n 100 e(i-1) = e(i) c f = 0.0e0 tst1 = 0.0e0 e(n) = 0.0e0 c do 240 l = 1, n j = 0 h = abs(d(l)) + abs(e(l)) if (tst1 .lt. h) tst1 = h c .......... look for small sub-diagonal element .......... do 110 m = l, n tst2 = tst1 + abs(e(m)) if (tst2 .eq. tst1) go to 120 c .......... e(n) is always zero, so there is no exit c through the bottom of the loop .......... 110 continue c 120 if (m .eq. l) go to 220 130 if (j .eq. 30) go to 1000 j = j + 1 c .......... form shift .......... l1 = l + 1 l2 = l1 + 1 g = d(l) p = (d(l1) - g) / (2.0e0 * e(l)) r = pythag(p,1.0e0) d(l) = e(l) / (p + sign(r,p)) d(l1) = e(l) * (p + sign(r,p)) dl1 = d(l1) h = g - d(l) if (l2 .gt. n) go to 145 c do 140 i = l2, n 140 d(i) = d(i) - h c 145 f = f + h c .......... ql transformation .......... p = d(m) c = 1.0e0 c2 = c el1 = e(l1) s = 0.0e0 mml = m - l c .......... for i=m-1 step -1 until l do -- .......... do 200 ii = 1, mml c3 = c2 c2 = c s2 = s i = m - ii g = c * e(i) h = c * p r = pythag(p,e(i)) e(i+1) = s * r s = e(i) / r c = p / r p = c * d(i) - s * g d(i+1) = h + s * (c * g + s * d(i)) c .......... form vector .......... do 180 k = 1, n h = z(k,i+1) z(k,i+1) = s * z(k,i) + c * h z(k,i) = c * z(k,i) - s * h 180 continue c 200 continue c p = -s * s2 * c3 * el1 * e(l) / dl1 e(l) = s * p d(l) = c * p tst2 = tst1 + abs(e(l)) if (tst2 .gt. tst1) go to 130 220 d(l) = d(l) + f 240 continue c .......... order eigenvalues and eigenvectors .......... do 300 ii = 2, n i = ii - 1 k = i p = d(i) c do 260 j = ii, n if (d(j) .ge. p) go to 260 k = j p = d(j) 260 continue c if (k .eq. i) go to 300 d(k) = d(i) d(i) = p c do 280 j = 1, n p = z(j,i) z(j,i) = z(j,k) z(j,k) = p 280 continue c 300 continue c go to 1001 c .......... set error -- no convergence to an c eigenvalue after 30 iterations .......... 1000 ierr = l 1001 return end real function pythag(a,b) real a,b c c finds sqrt(a**2+b**2) without overflow or destructive underflow c real p,r,s,t,u p = amax1(abs(a),abs(b)) if (p .eq. 0.0e0) go to 20 r = (amin1(abs(a),abs(b))/p)**2 10 continue t = 4.0e0 + r if (t .eq. 4.0e0) go to 20 s = r/t u = 1.0e0 + 2.0e0*s p = u*p r = (s/u)**2 * r go to 10 20 pythag = p return end real function epslon (x) real x c c estimate unit roundoff in quantities of size x. c real a,b,c,eps c c this program should function properly on all systems c satisfying the following two assumptions, c 1. the base used in representing floating point c numbers is not a power of three. c 2. the quantity a in statement 10 is represented to c the accuracy used in floating point variables c that are stored in memory. c the statement number 10 and the go to 10 are intended to c force optimizing compilers to generate code satisfying c assumption 2. c under these assumptions, it should be true that, c a is not exactly equal to four-thirds, c b has a zero for its last bit or digit, c c is not exactly equal to one, c eps measures the separation of 1.0 from c the next larger floating point number. c the developers of eispack would appreciate being informed c about any systems where these assumptions do not hold. c c this version dated 4/6/83. c a = 4.0e0/3.0e0 10 b = a - 1.0e0 c = b + b + b eps = abs(c-1.0e0) if (eps .eq. 0.0e0) go to 10 epslon = eps*abs(x) return end subroutine driver CVD\$G NOCONCUR c * real d(400),e(400),q(400,400) real alpha,alphap integer n,ldq common /prbdef/n,ldq,q,d,e c real z(400),beta(400),wv2(400),wm1(400,400),wv3(400) integer indx(400),indxp(400) common /wspace/z,beta,wv2,wm1,wv3,indx,indxp integer inodes(50),ndiml(50),ndimr(50) common/prfpms/ksect,kgran c mxtree = 50 off = 1 itql = 2 iupd = 3 ievu = 4 c c the test problem has been defined now comes the numerical c refinements that will avoid cancellation c do 20 i = 1,n z(i) = 0.0 20 continue nbp = n/2 inodes(1) = nbp + 1 ndiml(1) = nbp ndimr(1) = n - nbp klast = 1 jleft = 0 jright = 1 do 50 kk = 1,ksect-1 do 40 jj = 0,klast-1 jleft = jleft + 2 jright = jright + 2 klj = klast + jj ndiml(jleft) = ndiml(klj)/2 ndimr(jleft) = ndiml(klj) - ndiml(jleft) inodes(jleft) = inodes(klj) - ndimr(jleft) ndiml(jright) = ndimr(klj)/2 ndimr(jright) = ndimr(klj) - ndiml(jright) inodes(jright) = inodes(klj) + ndiml(jright) 40 continue klast = 2*klast 50 continue klast = 2**ksect j = klast klast = klast/2 do 90 jj = 1,klast j = j-1 ndl = ndiml(j) ir = inodes(j) il = ir - ndl nbp = ir-1 95 continue beta(nbp+1) = e(nbp+1) alpha = d(nbp) alphap = d(nbp+1) z(nbp) = 1.0e0 z(nbp+1) = 1.0e0 beta(nbp+2) = 1.0e0 if (sign(1.0e0,alpha)*sign(1.0e0,alphap) .ge. 0.0e0) * then if (alpha .lt. 0.0e0 .or. alphap .lt. 0.0e0) * then if(beta(nbp+1) .lt. 0.0e0) * then beta(nbp+1) = -beta(nbp+1) z(nbp) = 1.0e0 z(nbp+1) = -1.0e0 beta(nbp+2) = -1.0e0 endif else if (beta(nbp+1) .gt. 0.0e0) * then beta(nbp+1) = -beta(nbp+1) z(nbp) = 1.0e0 z(nbp+1) = -1.0e0 beta(nbp+2) = -1.0e0 endif endif endif d(nbp) = d(nbp) - beta(nbp+1) d(nbp+1) = d(nbp+1) - beta(nbp+1) if (nbp .eq. ir-1) * then nbp = il-1 if (nbp .gt. 0) go to 95 endif 90 continue c c invert the tree and spawn problems c klast = 2**ksect j = klast klast = klast/2 jtag = j + klast do 140 jj = 1,klast c c define the subproblems at the leaves c j = j-1 jtag = jtag - 1 c c place a problem on the queue here c the request is for call to tql on j-th leaf c to solve two tridiagonal eigenvalue problems c c myjob = itql icango = 0 nchks = 1 mchkin = j nparms = 3 indx(1) = ndiml(j) indx(2) = ndimr(j) indx(3) = inodes(j) call putq(jtag,myjob,icango,nchks,mchkin,nparms,indx) c ndl = indx(1) c ndr = indx(2) c ir = indx(3) c il = ir - ndl c c call tql2(ldq,ndl,d(il),e(il),q(il,il),ierr) c call tql2(ldq,ndr,d(ir),e(ir),q(ir,ir),ierr) c 140 continue klast = 2**ksect j = klast klast = klast/2 do 250 kk = 1,ksect do 240 jj = 1,klast j = j-1 c myjob = iupd icango = 1 if (klast .lt. 2**(ksect - 1) ) icango = 2 nchks = 1 mchkin = j/2 if(mchkin .eq. 0) nchks = 0 nparms = 3 indx(1) = ndiml(j) indx(2) = ndimr(j) indx(3) = inodes(j) call putq(j,myjob,icango,nchks,mchkin,nparms,indx) c 240 continue klast = klast/2 250 continue c return c c last card of driver c end subroutine sesupd(ldq,n,q,d,rho,z,x,dlamda,q2,delta,w,indxp, * indx,myj,ifail) CVD\$G NOCONCUR c*********************************************************************** c c c this subroutine will compute the updated eigensystem of a c of a symmetric matrix after modification by a rank one c symmetric matrix. c c a = qdq' + rho*z*z' c c it is assumed that the eigenvectors of the original matrix c are stored in q, and the eigenvalues are in d. c the algorithm consists of three stages... c c c the first stage constists of deflating the size of c the problem when there multiple eigenvalues or if there c zero of the vector q'z. for each such ocurrence the dimension c is reduced by one. c c the second stage consists of calculating the updated c eigenvalues of the reduced problem. this requires a call c to the zero finding routine evupd. c c the final stage consists of computing the updated eigenvectors c directly using the updated eigenvalue. c c c the algorithm requires o(n**2) operations to update the c eigenvectors, but n**3 + o(n**2) to update the eigenvectors. c c c input variables... c c n the dimension of the problem. q is n x n c c q an n x n matrix that contains the eigenvectors of c the original matrix on input and the updated c eigenvectors on output. c c d a vector of length n. the original eigenvalues are c contained in d on input. the updated eigenvectors c are contained on output. c c rho a scalar c c z a vector of length n. on input this vector c containes the updating vector. the contents of z c are destroyed during the updating process. c c x a working array of length n. c c dlamda a working array of length n c c q2 a working array of dimension n x n c c delta a working array of length n c c w a working array of length n. c c indx an integer array of length n. c c ifail this integer variable indicates failure of the c updating process with value 1, and success c with value 0. c c called subroutines... c c evupd a subroutine for calculating the updated eigenvalues. c c c*********************************************************************** real q(ldq,n),d(n),z(n),x(n),dlamda(n),q2(ldq,n),delta(n), 1 w(n) integer ldq,n,ifail integer indx(n),indxp(n) integer jndx(5) common/prfpms/ksect,kgran logical rhoge0,wait c c find out if this process has been reentered c and make appropriate branch c call enter(myj,label,k) if (label .gt. 0) go to 1111 c c c eps is machine precision c eps = epslon(1.0) zero = 0.0e0 ifail = 0 one = 1.0e0 two = 2.0e0 rhoge0 = (rho .ge. zero) c c compute q(transpose)*z c do 200 i = 1,n s = zero do 100 j = 1,n s = s + q(j,i)*z(j) q2(j,i) = q(j,i) 100 continue w(i) = s 200 continue c c normalize z so that norm(z) = 1 c t = enorm(n,w) do 300 j = 1,n z(j) = w(j)/t indx(j) = j 300 continue rho = rho*t*t if (rho .eq. zero) return c c order the eigenvalues c nm1 = n - 1 do 600 j = 1,nm1 t = d(1) s = z(1) inx = indx(1) k = n - j + 1 do 500 i = 2,k im1 = i - 1 if (d(i) .lt. t) go to 400 t = d(i) s = z(i) inx = indx(i) go to 500 400 continue d(im1) = d(i) d(i) = t z(im1) = z(i) z(i) = s indx(im1) = indx(i) indx(i) = inx 500 continue 600 continue c c calculate the spread in eigenvalues c evsprd = abs(d(1)-d(n)) c c if there multiple eigenvalues then the problem deflates. c here the number of equal eigenvalues are found. c then an elementary reflector is computed to rotate the c corresponding eigensubspace so that certain components of c z are zero in this new basis. c if (rhoge0) *then j = 1 dlam = d(1) k = 0 k2 = n + 1 t = z(1)**2 icount = 1 700 continue jm1 = j j = j + 1 if (j .gt. n) go to 800 if(abs(d(j)-dlam) .le. eps*evsprd) * then t = t + z(j)**2 icount = icount + 1 k2 = k2 - 1 indxp(k2) = jm1 else dmax = amax1(abs(d(1)),abs(d(n))) if (abs(rho)*sqrt(t) .gt. n*dmax*eps) * then k = k + 1 dlamda(k) = d(jm1) indxp(k) = jm1 if (icount .gt. 1) * call rotate(icount,indx,jm1,n,ldq,q,z,rhoge0) w(k) = t else k2 = k2 - 1 indxp(k2) = jm1 endif dlam = d(j) icount = 1 t = z(j)**2 endif go to 700 800 continue dmax = amax1(abs(d(1)),abs(d(n))) if (abs(rho)*sqrt(t) .gt. n*dmax*eps) * then k = k + 1 w(k) = t dlamda(k) = d(jm1) indxp(k) = jm1 if (icount .gt. 1) * call rotate(icount,indx,jm1,n,ldq,q,z,rhoge0) else k2 = k2 - 1 indxp(k2) = jm1 endif else jm1 = n dlam = d(n) k = n + 1 k2 = 0 t = z(n)**2 icount = 1 701 continue j = jm1 jm1 = j - 1 if (jm1 .lt. 1) go to 801 if(abs(d(jm1)-dlam) .le. eps*evsprd) * then t = t + z(jm1)**2 icount = icount + 1 k2 = k2 + 1 indxp(n-k2+1) = j else dmax = amax1(abs(d(1)),abs(d(n))) if (abs(rho)*sqrt(t) .gt. n*dmax*eps) * then k = k - 1 w(n-k+1) = t dlamda(n-k+1) = -d(j) indxp(n-k+1) = j jp = jm1 + icount if (icount .gt. 1) * call rotate(icount,indx,jp,n,ldq,q,z,rhoge0) else k2 = k2 + 1 indxp(n-k2+1) = j endif dlam = d(jm1) icount = 1 t = z(jm1)**2 endif go to 701 801 continue dmax = amax1(abs(d(1)),abs(d(n))) if (abs(rho)*sqrt(t) .gt. n*dmax*eps) * then k = k - 1 w(n-k+1) = t dlamda(n-k+1) = -d(j) indxp(n-k+1) = j jp = jm1 + icount if (icount .gt. 1) * call rotate(icount,indx,jp,n,ldq,q,z,rhoge0) else k2 = k2 + 1 indxp(n-k2+1) = j endif k = n - k + 1 rho = -rho endif c c c compute the updated eigenvalues of the deflated problem c c nmax = 2**ksect - 1 nproc = 1 ksize = k node = myj idepth = 1 2640 continue if (ksize .gt. 2*kgran .and. node .lt. nmax) then nproc = nproc*2 node = node*2 ksize = k/nproc + 1 idepth = idepth + 1 go to 2640 endif kstart = k - ksize + 1 kstop = k jleft = myj c c prepare to spawn processes c if ( kstart .gt. ksize ) call prtspn(myj) kproc = 0 c do 2660 kk = 2,idepth jleft = jleft*2 do 2670 jj = jleft,jleft + 2**(kk-1) - 1 if (kstart .gt. ksize) * then myjob = 4 icango = 0 nparms = 5 jndx(1) = k jndx(2) = kstart jndx(3) = kstop jndx(4) = rhoge0 jndx(5) = myj call spawn(myj,myjob,nparms,jndx) c call evdrv(k,kstart,kstop,ldq,n,q,d,rho,z,x, c * dlamda,q2,delta,w,indxp,indx,rhoge0,ifail) kproc = kproc + 1 kstop = kstart - 1 kstart = kstop - ksize + 1 endif 2670 continue 2660 continue kstart = 1 call evdrv(k,kstart,kstop,ldq,n,q,d,rho,z,x, * dlamda,q2,delta,w,indxp,indx,rhoge0,ifail) c c if any processes were spawned then c place a reentry point here and resume when all spawned c processes have checked in c if( kproc .gt. 0 .and. wait(myj,1,1,k)) return 1111 continue c c c store the updated eigenvectors back into q c do 2780 j = k+1,n jp = indxp(j) jjpp = indx(jp) do 2740 i = 1,n q2(i,jp) = q(i,jjpp) 2740 continue 2780 continue do 2900 j = 1,n do 2800 i = 1,n q(i,j) = q2(i,j) 2800 continue 2900 continue c ifail = n-k c return c c last card of sesupd c end subroutine evdrv(k,kstart,kstop,ldq,n,q,d,rho,z,x,dlamda,q2, * delta,w,indxp,indx,rhoge0,ifail) CVD\$G NOCONCUR c c real q(ldq,n),d(n),z(n),x(n),dlamda(n),q2(ldq,n),delta(n), * w(n),rho integer k,kstart,kstop,ldq,n,ifail integer indx(n),indxp(n) logical rhoge0 c c local variables c real t,s,dlam c do 2700 j = kstart,kstop i = j if (.not. rhoge0) i = k - i + 1 call evupd(k,i,dlamda,w,delta,rho,dlam,ifail) c c if the zero finder failed the computation is terminated c if (ifail .eq. 1) return c if (rhoge0) * then jp = indxp(j) else dlam = -dlam jp = indxp(k-j+1) endif c d(jp) = dlam c c compute the updated eigenvectors c do 2200 i1 = 1,n x(i1) = 0.0e0 2200 continue do 2500 jj = 1,k j3 = jj if (.not. rhoge0) j3 = k - jj + 1 jp = indxp(j3) jjpp = indx(jp) s = z(jp)/delta(j3) do 2300 i1 = 1,n x(i1) = x(i1) + q(i1,jjpp)*s 2300 continue 2500 continue t = enorm(n,x) j3 = j if (.not. rhoge0) j3 = k - j + 1 jp = indxp(j3) do 2600 i1 = 1,n q2(i1,jp) = x(i1)/t 2600 continue 2700 continue c return c c last card of evdrv c end subroutine work(id) CVD\$G NOCONCUR c****************************************************************************** c c this subroutine is a template for scheduling parallel processing c the user is responsible for declaring ALL shared variables in common c area common/prbdef/ and for declaring scratch space of sufficient c size for any schedulable process associated with the task served by c work. c c****************************************************************************** c * c user specifies probdef which contains all shared variables * c * real d(400),e(400),q(400,400) integer n,ldq common /prbdef/n,ldq,q,d,e c real z(400),beta(400),wv2(400),wm1(400,400),wv3(400) integer indx(400),indxp(400) common /wspace/z,beta,wv2,wm1,wv3,indx,indxp common /histry/ihist(100,8) c c * c****************************************************************** c c declare local storage c real ww1(400),ww2(400) integer gtprb integer jndx(50) c c********************************************** c * c user specifies all scratch arrays here * c * real wv4(400) , wv1(400) logical rhoge0 c * c********************************************** c c worker with id = 1 is designated master c others are slaves c if (id .eq. 1) then c c********************************************** c * c the user must supply driver to * c define schedulable processes and * c data dependencies. * call driver c * c********************************************** call start endif c kd = 1 c 10 continue myjob = gtprb(jobtag) ihist(kd,id) = jobtag kd = kd + 1 call gtprms(jobtag,jndx) c********************************************** c * c the user must make myjob specified * c in driver correspond to the label * c of the appropriate subroutine call * c this should be an integer .ge. 2 * c since label 100 is reserved for * c return code. * c * c go to (100,200,300,400),myjob c * c********************************************** c c 100 continue c return c 200 continue c c********************************************** c * c the user must make myjob request * c correspond to this subroutine call * c the parameters which define the * c call are indices that select areas * c from the approriate arrays in /probdef/ * c * c request is itql c ndl = jndx(1) ndr = jndx(2) ir = jndx(3) il = ir - ndl c call tql2(ldq,ndl,d(il),e(il),q(il,il),ierr) call tql2(ldq,ndr,d(ir),e(ir),q(ir,ir),ierr) c c * c********************************************** call chekin(jobtag) go to 10 c 300 continue c c********************************************** c c request is iupd c ndl1 = jndx(1) ndr1 = jndx(2) ir1 = jndx(3) il1 = ir1 - ndl1 ndim1 = ndl1+ndr1 c c n q d beta z are inputs c z(ir1-1) = 1.0e0 z(ir1) = 1.0e0*beta(ir1+1) do 222 jjj = il1,ir1-2 z(jjj) = 0.0e0 222 continue do 333 jjj = ir1+1,ir1+ndr1-1 z(jjj) = 0.0e0 333 continue myj = jobtag c call sesupd(ldq,ndim1,q(il1,il1),d(il1),beta(ir1),z(il1), * wv1(il1),wv2(il1),wm1(il1,il1),wv4(il1),wv3(il1), * indxp(il1),indx(il1),myj,ifail) c c********************************************** c call chekin(jobtag) go to 10 c 400 continue c c********************************************** c c request was for ievu c k = jndx(1) kstart = jndx(2) kstop = jndx(3) rhoge0 = jndx(4) mypar = jndx(5) c call gtprms(mypar,jndx) ndl = jndx(1) ndr = jndx(2) ir = jndx(3) il = ir - ndl ndim = ndl+ndr c call evdrv(k,kstart,kstop,ldq,ndim,q(il,il),d(il),beta(ir), * z(il),ww2(il),wv2(il),wm1(il,il),ww1(il),wv3(il), * indxp(il),indx(il),rhoge0,ifail) c********************************************** c call chekin(jobtag) go to 10 c 500 continue c c********************************************** c additional sub call placed here c extend this pattern further if needed. c********************************************** c call chekin(jobtag) go to 10 c end