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