program opsol
         USE linear_operators
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.
c	call rand_gen(y)
         A = rand(A)
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.

         b = rand(b)
c	write(6,100)((b(i,j),j=1,1),i=1,n)
! Compute the solution matrix of Ax=b.
         t1=etime(dum)
         x=A.ix.b
         t2=etime(dum)
         res=b-(A.x.x)
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. 
c	res = b - matmul(A,x)
c	write(6,200)((res(i,j),j=1,1),i=1,n)
c200 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