c Project a set of np vectors from n+1 coords summing to 0 to n coords. c Uses standard input/output. c The print statement in line 33 is set up to print four coordinates c per line. double precision a(30,30),bef(100),aft(100,30) double precision x,sum,scale nmax=30 npmax=100 c read n (the real dimension) and np, the number of points read(*,*)n,np if(n.le.nmax.and.np.le.npmax)goto 570 write(*,*)" n or np too big" stop 570 continue c construct the projection matrix, a, n+1 x n do 10 i=1,n+1 do 10 j=1,n 10 a(i,j)=0.0d0 do 11 j=1,n x=dfloat(j*(j+1)) x=1.0d0/dsqrt(x) do 12 i=1,j 12 a(i,j)=x 11 a(j+1,j)=-dfloat(j)*x write(*,*)" a" do 13 i=1,n+1 13 write(*,100)(a(i,j),j=1,n) write(*,*)"before" c read np vectors each with np+1 coordinates do 14 ip=1,np read(*,*)(bef(i),i=1,n+1) write(*,100)(bef(i),i=1,n+1) 100 format(4d28.18) c check sum is 0 sum=0.0d0 do 1 i=1,n+1 1 sum=sum+bef(i) if(dabs(sum).lt.0.00001d0)goto 2 write(*,*)"Error: sum not 0 ",sum stop 2 continue do 15 i=1,n sum=0.0d0 do 17 j=1,n+1 17 sum=sum+bef(j)*a(j,i) aft(ip,i)=sum 15 continue c normalize to have length 1 sum=0.0d0 do 19 i=1,n 19 sum=sum+aft(ip,i)**2 scale=1.0d0/dsqrt(sum) do 20 i=1,n 20 aft(ip,i)=scale*aft(ip,i) 14 continue write(*,*)"after" do 16 ip=1,np write(*,100)(aft(ip,i),i=1,n) 16 continue stop end