c**************************************************************

	SUBROUTINE least_sq(a,b,m,n,x,eps,p)

c Routine PUBLIC de mod_numerique	
c Résolution d'un système linéaire surdéterminé par moindres carrés

c entrées:
c	a(m,n) : coefficients du système
c	b(m) : second membbre
c	m : nombre de lignes (observations)
c	n : nombre de colonnes (inconnues)
c	p : poids

c sorties:
c	x : solution
c	eps : précison

c Auteur: P. Morel, Département Cassiopée, O.C.A.

c----------------------------------------------------------------

	USE mod_kind
	
	IMPLICIT NONE
	
	REAL (kind=dp), INTENT(in), DIMENSION(m,n) :: a
	REAL (kind=dp), INTENT(in), OPTIONAL, DIMENSION(m) :: p
	REAL (kind=dp), INTENT(in), DIMENSION(m) :: b
	INTEGER, INTENT(in) :: m,n
	REAL (kind=dp), INTENT(out), DIMENSION(n) :: eps, x
	
	REAL (kind=dp), DIMENSION(n,n) :: an
	REAL (kind=dp), DIMENSION(n,m) :: at
	REAL (kind=dp), DIMENSION(m,n) :: pa
	REAL (kind=dp), DIMENSION(m) :: pb, res
	REAL (kind=dp), DIMENSION(n) :: bi
	REAL (kind=dp) :: sig
	
	INTEGER :: i
	
	LOGICAL :: ok
	
c---------------------------------------------------------------

2000	FORMAT(8es10.3)

c P A  et P B
	IF(PRESENT(p))THEN
	 pb=p*b
	 DO i=1,m
	  pa(i,:)=p(i)*a(i,:)
	 ENDDO
	ELSE
	 pb=b ; pa=a
	ENDIF
	
c système linéaire At P A = At P B
	at=TRANSPOSE(a) ; an=MATMUL(at,pa) ; CALL matinv(an,n,ok)	
	bi=MATMUL(at,pb) ; x=MATMUL(an,bi)

c residus A X - B
	res=(MATMUL(a,x)-b)**2

	IF(PRESENT(p))THEN	
c sig2=SUM P * res2 / m-n
	 sig=SQRT(SUM(p*res)/REAL(m-n,dp))
	 
	ELSE
c sig2=SUM res2 / m-n
	 sig=SQRT(SUM(res)/REAL(m-n,dp))
	ENDIF

c précision		
	DO i=1,n
	 eps(i)=sig*SQRT(an(i,i))
	ENDDO
		
	RETURN
	
	END SUBROUTINE least_sq