double precision dd(400),ee(400),qq(400,400) double precision esave(400),dsave(400) double precision s,s2,t,t2,tnorm,tnorm2,res,res2 integer icase c double precision d(400),e(400),q(400,400) integer n,ldq common /prbdef/n,ldq,q,d,e common/prfpms/ksect,kgran c c ldq = 400 c read(5,*) nproc read(5,*) ksect read(5,*) kgran write(6,*) ' input ' write(6,*) nproc,ksect,kgran do 9999 n = 50,100,50 write(6,*)'==============================' write(6,*)' n = ',n do 9998 icase = 1,2 write(6,*)'++++++++++++++++++++++++++++++' if( icase .eq. 1 ) write(6,*)' twos on diagonal' if( icase .eq. 2 ) write(6,*)' random numbers on diagonal' nsect = 2**ksect nd2 = n/2 go to (112,113), icase 112 do 13 i = 1,n d(i) = 2.0 e(i) = 1 13 continue go to 444 113 do 14 i = 1,n d(i) = rand(foo) e(i) = rand(foo) 14 continue 444 do 10 j = 1,n do 5 i = 1,n q(i,j) = 0.0 qq(i,j) = 0.0 5 continue q(j,j) = 1.0 qq(j,j) = 1.0 dd(j) = d(j) ee(j) = e(j) dsave(j) = d(j) esave(j) = e(j) 10 continue c if (n .eq. 50) go to 9999 t1 = second(gtime) c call tql2(ldq,n,dd,ee,qq,ierr) t2t = second(gtime) - t1 write(6,*) ' time for tql ',t2t c c the test problem has been defined now comes the numerical c refinements that will avoid cancellation c t1 = second(gtime) c call libopn(nproc) c t2 = second(gtime) - t1 write(6,*) ' time for sesupd ',t2 if( ifail .gt. 1 ) write(6,*)' deflate from sesupd',ifail write(6,*)' ratio of tql2/new',t2t/t2 tnorm = 0.0d0 tnorm2 = 0.0d0 res = 0.0d0 res2 = 0.0d0 do 530 j = 1,n e(1) = dsave(1)*q(1,j) + \$ esave(2)*q(2,j) - d(j)*q(1,j) ee(1) = dsave(1)*qq(1,j) + \$ esave(2)*qq(2,j) - dd(j)*qq(1,j) do 400 i = 2,n-1 e(i) = esave(i)*q(i-1,j) + dsave(i)*q(i,j) + \$ esave(i+1)*q(i+1,j) - d(j)*q(i,j) ee(i) = esave(i)*qq(i-1,j) + dsave(i)*qq(i,j) + \$ esave(i+1)*qq(i+1,j) - dd(j)*qq(i,j) 400 continue e(n) = esave(n)*q(n-1,j) + dsave(n)*q(n,j) \$ - d(j)*q(n,j) ee(i) = esave(n)*qq(n-1,j) + dsave(n)*qq(n,j) \$ - dd(j)*qq(n,j) t = enorm(n,e) t2 = enorm(n,ee) res = max(res,t) res2 = max(res2,t2) if (t .gt. .01) then write(6,*)' j ',j, ' d ',d(j),' er ',t,' dd ',dd(j),' eer ',t2 write(6,*) ' ev ',q(1,j),' eev ',qq(1,j) endif do 520 i = 1,n t = 0.0 s = 0.0 t2 = 0.0 s2 = 0.0 do 510 k = 1,n s2 = s2 + qq(k,i)*qq(k,j) s = s + q(k,i)*q(k,j) 510 continue if (i .eq. j) s2 = s2 - 1.0 if (i .eq. j) s = s - 1.0 t2 = abs(s2) t = abs(s) tnorm2 = max(tnorm2,t2) tnorm = max(tnorm,t) c if (t .gt. .000000001) write(6,*) ' ij ',i,j,' tnorm ',tnorm 520 continue 530 continue write(6,*)' the residual for the tql values and vectors',res2 write(6,*)' the residual for the updated values and vectors',res write(6,*)' the tql norm of q*q sup t is',tnorm2 write(6,*)' the upd norm of q*q sup t is',tnorm 9998 continue 9999 continue stop end