! --------------------------------------------------------------------------- !
!     Test program for LINEAR_SOLVE [real and quadruple version]              !
!       A:=sqrt[2/(n+1)] sin[(i*j*pi)/((N+1)]                                 !
!             -1                                                              !
!         --> A = A                                                           !
!             The solution of Ax=b is obtained by A*b                         !
! --------------------------------------------------------------------------- !
program linear_solve_test
 use GECP
 implicit none

 integer                                    :: n
 real(kind(1Q0)),dimension(:,:),allocatable :: A
 real(kind(1Q0)),dimension(:),allocatable   :: x1,x2,b,sol
 integer                                    :: i,j
 real(kind(1Q0))                            :: Pi,s

 n=100
 allocate(A(n,n),x1(n),x2(n),b(n),sol(n))

 Pi=4*atan(1.0Q0)
 A=0
 do j=1,n
   do i=1,n
     A(i,j)=sqrt(2.0Q0/(n+1))*sin((i*j*Pi)/(n+1))
   end do
 end do

 call random_number(b)            
 sol=matmul(A,b)

 x1=linear_solve(A,b)                     ! solve linear equation
! x2=linear_solve(A,b,sv=s)          

 write(6,*) 
 write(6,*) 'dimension=',N
 write(6,*) '               maximum of the solution=',maxval(x1)
 write(6,*) '        relative error of linear solve=',&
                    maxval(abs(x1-sol))/maxval(abs(sol))
 write(6,*) ' singular value=',s

 deallocate(A,x1,x2,b,sol)

end program linear_solve_test