program sol90
         use lin_sol_gen_int
         use rand_gen_int
         use error_option_packet

c	implicit none

! This is Example 1 for LIN_SOL_GEN.

         integer, parameter :: n=1024
         double precision, parameter :: one=1d0
         double precision err
         double precision A,b,x,res
         dimension A(n,n), b(n,1), x(n,1), res(n,1) real(kind(1e0)) y(n**2)
c	real*8 s
! Generate a random matrix.
         call rand_gen(y)
         A = reshape(y,(/n,n/))
c	write(6,100)((A(i,j),j=1,n),i=1,n)
c100 format(8e10.3)
         ops=2*n*n*n/3.0
! Generate random right-hand sides.
         call rand_gen(y)
         b = reshape(y,(/n,1/))
c	write(6,100)((b(i,j),j=1,1),i=1,n)
! Compute the solution matrix of Ax=b.
         t1=etime(dum)
         call lin_sol_gen(A, b, x)
c	do 15 i=1,n
c	do 15 j=1,n
c	s=0.0d0
c	do 5 k=1,n
c	s=s+a(i,k)*x(k,j)
c 5	continue
c	res(i,j)=s-b(i,j)
c 15	continue
         t2=etime(dum)
c	write(6,100)((x(i,j),j=1,1),i=1,n)
         write(*,*)'n=',n,t2-t1,'sec ',ops/(t2-t1)*1.e-6,'Mflops' 
! Check the results for small residuals. 
         res = b - matmul(A,x)
         write(6,200)((res(i,j),j=1,1),i=1,n)
  200 format(4d15.8)
         err = maxval(abs(res))/(sum(abs(A))+sum(b)) 
         if (err <= sqrt(epsilon(one))) then
           write (*,*) 'Example 1 for LIN_SOL_GEN is correct.' 
         end if

         end
         
         function etime( )
         etime=timef()*0.001
         return
         end