program test_lstsq
implicit none

double precision, parameter :: PI = 3.141592, C1=-4.87e2, C2=5.3411e3, C3=1.61153
double precision, allocatable, dimension(:,:) :: a, b
double precision, allocatable, dimension(:) :: t, work
integer, allocatable, dimension(:) :: iwork
real, dimension(10) :: singular
integer :: i, j, N, rank, lwork, info

allocate(iwork(3*3*1 + 11*3))

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
    
    allocate(work(1))
    call dgelsd(N, 3, 1, a, N, b, N, singular, -1.0, rank, work, -1, iwork, info)
    lwork = work(1)
    
    deallocate(work)
    allocate(work(lwork))
    
    singular = 0
    rank = 0
    work = 0
    iwork = 0
    
    call dgelsd(N, 3, 1, a, N, b, N, singular, -1.0, rank, work, lwork, iwork, info)
    deallocate(work)
    
    write(*,'(A, I6, A5, 3E15.3)') 'N=', N, 'X=', b(:3,1) 
    
    !write(*,*) N
    !write(*,*) 'answers', b(:3,1)
    !write(*,*) 'singular', singular 
    !write(*,*) lwork, rank, info
    !write(*,*) 
    
    deallocate(a)
    deallocate(b)
    deallocate(t)
    N = N*2
end do

deallocate(iwork)

end program
