program test_lstsq
use f95_lapack, only: la_gelsy, la_gelsd
implicit none

double precision, parameter :: PI = 4*tan(1.0), C1=-4.87e2, C2=5.3411e3, C3=1.61153
double precision, allocatable, dimension(:,:) :: a, b
double precision, allocatable, dimension(:) :: t
integer :: i, j, N, rank, info

write(*,*) C1, C2, C3
write(*,*)

N = 100
do i = 1, 8
    allocate(a(N,3))
    allocate(b(N,1))
    allocate(t(N))
    
    do j = 1, N
        t(j) = j*2*PI/N
    end do
    
    b(:,1) = C1*sin(t) + C2*cos(t) + C3
    a(:,1) = sin(t)
    a(:,2) = cos(t)
    a(:,3) = 1
    
    rank = 0
    info = 0
    
    call la_gelsd(a, b, rank=rank, info=info)
    
    write(*,'(A, I6, A5, 3E15.3)') 'N=', N, 'X=', b(:3,1) 
    
    deallocate(a)
    deallocate(b)
    deallocate(t)
    N = N*2
end do


end program
