program vande
c	***************************************************
c	Inverts a Vandermonde matrix explicitly
c	n    :  dimension of matrix to be inverted
c	iflg :  0=matrix inversion only; 1=equation solving
c	m    :  number of invariants found using the Newton identities;
c               the remaining n-m are found using the conjugate
c               Newton identities: (0.le.m.le.n)
c	        Recommended value for lamda(i)=i: m=2n/3
c      Redirect output to file: Example  vande.e >> vande.out
c	***************************************************
	integer i,j,k,n,m
	real*8 a(100,100),inv(100,100),an(100,100),trn(100),
     &         tr(100),f(100),x(100),b(100),lamda(100),den(100),
     &         trn1(100)
	open(unit=5,file='vande.dat',status='old')
	open(unit=6,file='vande.out',status='new')
	read(5,*) n,iflg,m
	read(5,*) (lamda(i),i=1,n)
	if(iflg.eq.1)then
	read(5,*) (f(i),i=1,n)
	endif
	if(m.lt.0.or.m.gt.n)then
	 write(6,*) 'm has to lie between 0 and n '
	 stop
	endif
	do 998 i=1,n
	 if(lamda(i).eq.0.0d0)then
	  m=n-1
	  go to 999
	 endif
998	continue
999	continue
	do 1 i=1,n
	do 1 j=1,n
	 a(i,j)=lamda(j)**(i-1)
1	continue
	do 3 i=1,n
	 den(i)=1.0d0
	do 2 j=1,n
	 if(i.ne.j)then
	 den(i)=den(i)*(lamda(i)-lamda(j))
	 endif
2	continue
	if(den(i).eq.0.0d0)then
	write(6,*) 'Matrix not invertible'
        stop
        endif
3	continue
	call zeror(inv,100*100)
	do 100 i=1,n
	 do 7 j=1,n
	  tr(j)=lamda(j)
	  trn1(j)=lamda(j)
7	 continue
	tr(i)=0.0d0
	if(m.gt.0)then
	inv(i,n)=1.0d0
	endif
	if(m.gt.1)then
	do 8 j=1,n
	inv(i,n-1)=inv(i,n-1)-tr(j)
8	continue
	endif
	do 24 j=2,m-1
	 do 4 k=1,n
	  trn(k)=trn1(k)
4	 continue
	  do 17 k=1,n
	   trn(k)=trn(k)+inv(i,n-j+1)
17	  continue
	 do 18 k=1,n
	  trn1(k)=trn(k)*tr(k)
18	 continue
	do 22 k=1,n
	inv(i,n-j)=inv(i,n-j)-trn1(k)
22	continue
	inv(i,n-j)=inv(i,n-j)/dfloat(j)
24	continue
	if(m.lt.n)then
	inv(i,1)=dfloat((-1)**(n-1))
	do 34 k=1,n
	 if(k.ne.i)then
	 inv(i,1)=inv(i,1)*lamda(k)
	 endif
34	continue 
	endif
	 if(m.lt.n-1)then
        do 35 k=1,n
         if(k.ne.i)then
         inv(i,2)=inv(i,2)-inv(i,1)/lamda(k)
         endif
35      continue
	 do 70 j=1,n
	  tr(j)=1.0/lamda(j)
	  trn1(j)=inv(i,1)/lamda(j)
70	 continue
	tr(i)=0.0d0
	do 240 j=2,n-m-1
	 do 40 k=1,n
	  trn(k)=trn1(k)
40	 continue
	  do 170 k=1,n
	   trn(k)=trn(k)+inv(i,j)
170	  continue
	 do 180 k=1,n
	  trn1(k)=trn(k)*tr(k)
180	 continue
	do 220 k=1,n
	inv(i,j+1)=inv(i,j+1)-trn1(k)
220	continue
	inv(i,j+1)=inv(i,j+1)/dfloat(j)
240	continue
	endif
100	continue
c	********************
c	compute the inverse
c	********************
	do 90 j=1,n
	do 90 k=1,n
	inv(j,k)=inv(j,k)/den(j)
90	continue
	if(iflg.eq.0)then
	write(6,*) 'Inverse of the matrix'
	write(6,*)
	do 120 j=1,n
	write(6,110) (inv(j,k),k=1,n)
110	format(10e13.5)
120	continue
c	*********************************
c	check the accuracy of the inverse
c	*********************************
	write(6,*) 
	write(6,*) 'Matrix*Inverse'
	write(6,*)
	call matmult(n,n,n,a,inv,an,100,100,100,100,100,100)
	do 112 j=1,n
	write(6,111) (an(j,k),k=1,n)
111	format(10e13.5)
112	continue
	endif
c	**********************************
c	Find the solution x=a^{-1}f
c	**********************************
	if(iflg.eq.1)then
	do 135 i=1,n
	x(i)=0.0d0
	do 135 k=1,n
	 x(i)=x(i)+inv(i,k)*f(k)
135	continue
	do 150 i=1,n
	b(i)=0.0d0
        do 150 j=1,n
        b(i)=b(i)+a(i,j)*x(j)
150     continue
	write(6,*) 
	write(6,*) 'Soln'
	write(6,*) 
	write(6,*) (x(i),i=1,n)
138     format(6e13.5)
	write(6,*) 
	write(6,*) 'Matrix*Soln'
	write(6,*)
        write(6,*) (b(i),i=1,n)
140     format(6e13.5)
	endif
	rewind(5)
	stop
	end
c       *********************************************************
c       *********************************************************
        subroutine matmult(l,m,n,a,b,c,ndim1,ndim2,ndim3,ndim4,
     1  ndim5,ndim6)
c       *********************************************************
c       subroutine multiplying 2 matrices a and b  to
c       give product matrix c  :  c(l,n)=a(l,m)*b(m,n)
c       *********************************************************
        dimension a(ndim1,ndim2),b(ndim3,ndim4),c(ndim5,ndim6)
        real*8 a,b,c,sum
        integer i,j,k,l,m,n,ndim1,ndim2,ndim3,ndim4,ndim5,ndim6
        do 10 i=1,l
        do 10 k=1,n
10      c(i,k)=0.0d0
        do 20 i=1,l
        do 20 k=1,n
        sum=0.0d0
        do 20 j=1,m
        sum=sum+a(i,j)*b(j,k)
20      c(i,k)=sum
        return
        end
c       *********************************************************
c       *********************************************************
        subroutine zeror (vector,isize)
c       *********************************************************
        dimension vector(isize)
        integer   isize
        real*8    vector
c
c       Local variables
c
        integer   i

        do 10 i = 1 , isize
          vector(i) = 0.0d0
10      continue

        return
        end