program lsq
implicit none
!
!This program calculates least square solutions for a model with 2 effects.
!All the storage is dense, and solutions are obtained iteratively by
!Gauss-Seidel.
! This program is easily upgradable to any number of effects
!
real, allocatable:: xx(:,:),xy(:),sol(:) !storage for the equations
integer, allocatable:: indata(:) !storage for one line of effects
!integer,parameter:: neff=2,nlev(2)=(/2,3/) !number of effects and levels
real :: y ! observation value
integer :: neq,io,i,j ! number of equations and io-status
integer :: neff,nlev(50)
integer,allocatable::address(:)
character (40) :: parfile,datafile
print*,'Parameter file?'
read(*,'(a)') parfile
open(99,file=parfile)
nlev=0
read(99,'(a)') datafile
print*,'Data file: ',datafile
read(99,*) neff
do
read(99,*,iostat=io) (nlev(i),i=1,neff)
if(io /= 0) exit
enddo
print*,'# effects =',neff
!
neq=sum(nlev)
print*,'# equations =',neq
allocate (xx(neq,neq), xy(neq), sol(neq),indata(neff),address(neff))
xx=0; xy=0; sol=0
!
open(1,file=datafile)
!
do
read(1,*,iostat=io)indata,y
if (io.ne.0) exit
call find_addresses
do i=1,neff
do j=1,neff
xx(address(i),address(j))=xx(address(i),address(j))+1
enddo
xy(address(i))=xy(address(i))+y
enddo
enddo
!
print*,'left hand side'
do i=1,neq
print '(100f5.1)',xx(i,:)
enddo
!
print '( '' right hand side:'' ,100f6.1)',xy
!
call solve_dense_gs(neq,xx,xy,sol) !solution by Gauss-Seidel
print '( '' solution:'' ,100f7.3)',sol
contains
subroutine find_addresses
integer :: i
do i=1,neff
address(i)=sum(nlev(1:i-1))+indata(i)
enddo
end subroutine
end program lsq
subroutine solve_dense_gs(n,lhs,rhs,sol)
! finds sol in the system of linear equations: lhs*sol=rhs
! the solution is iterative by Gauss-Seidel
integer :: n
real :: lhs(n,n),rhs(n),sol(n),eps
integer :: round
!
round=0
do
eps=0; round=round+1
do i=1,n
solnew=sol(i)+(rhs(i)-sum(lhs(i,:)*sol))/lhs(i,i)
eps=eps+ (sol(i)-solnew)**2
sol(i)=solnew
end do
if (eps.lt. 1e-10) exit
end do
print*,'solutions computed in ',round,' rounds of iteration'
end subroutine