subroutine Least_SQ_Fit(xo,yo,kj,mmx,mmy,mmxy,stxz,styz,jac,MATAX
     &,MATAY,MATAXY,MATAXZ,MATAYZ)
	use msimsl
	implicit none
	include  'eldata.h'
	include  'cdata.h'
	include  'sdata.h'
	include  'iofile.h'
	include  'pointer.h'
	include  'upointer.h'
	include  'comblk.h'

	real*8	mmx(16),mmy(16),mmxy(16),stxz(16),styz(16)
	real*8	xo(16),yo(16),ap(8,8),iap(8,8),jac(16)
	real*8   BMATX(8),BMATY(8),BMATXY(8),BMATXZ(8),BMATYZ(8)
	real*8   MATAX(8,numnp),MATAY(8,numnp),MATAXY(8,numnp)
     &        ,MATAXZ(8,numnp),MATAYZ(8,numnp),matp(128)
	integer i,j,ij,ji,kj,ii,ik,kk,l,is,jj
	
	do i=1,8
	do j=1,8
	ap(i,j)=0.0d0
	enddo
	enddo

	do i=1,8
	BMATX(i)=0.0d0
	BMATY(i)=0.0d0
	BMATXY(i)=0.0d0
	BMATXZ(i)=0.0d0
	BMATYZ(i)=0.0d0
	enddo

	do ii=1,16
	ap(1,1)=ap(1,1)+1.0d0
	ap(1,2)=ap(1,2)+xo(ii)
	ap(1,3)=ap(1,3)+yo(ii)
	ap(1,4)=ap(1,4)+(xo(ii)*xo(ii))
	ap(1,5)=ap(1,5)+(xo(ii)*yo(ii))
	ap(1,6)=ap(1,6)+(yo(ii)*yo(ii))
	ap(1,7)=ap(1,7)+(xo(ii)*xo(ii)*yo(ii))
	ap(1,8)=ap(1,8)+(xo(ii)*yo(ii)*yo(ii))
	ap(2,2)=ap(2,2)+(xo(ii)*xo(ii))
	ap(2,3)=ap(2,3)+(xo(ii)*yo(ii))
	ap(2,4)=ap(2,4)+(xo(ii)*xo(ii)*xo(ii))
	ap(2,5)=ap(2,5)+(xo(ii)*xo(ii)*yo(ii))
	ap(2,6)=ap(2,6)+(xo(ii)*yo(ii)*yo(ii))
	ap(2,7)=ap(2,7)+(xo(ii)*xo(ii)*xo(ii)*yo(ii))
	ap(2,8)=ap(2,8)+(xo(ii)*xo(ii)*yo(ii)*yo(ii))
	ap(3,3)=ap(3,3)+(yo(ii)*yo(ii))
	ap(3,4)=ap(3,4)+(xo(ii)*xo(ii)*yo(ii))
	ap(3,5)=ap(3,5)+(xo(ii)*yo(ii)*yo(ii))
	ap(3,6)=ap(3,6)+(yo(ii)*yo(ii)*yo(ii))
	ap(3,7)=ap(3,7)+(xo(ii)*xo(ii)*yo(ii)*yo(ii))
	ap(3,8)=ap(3,8)+(xo(ii)*yo(ii)*yo(ii)*yo(ii))
	ap(4,4)=ap(4,4)+(xo(ii)*xo(ii)*xo(ii)*xo(ii))
	ap(4,5)=ap(4,5)+(xo(ii)*xo(ii)*xo(ii)*yo(ii))
	ap(4,6)=ap(4,6)+(xo(ii)*xo(ii)*yo(ii)*yo(ii))
	ap(4,7)=ap(4,7)+((xo(ii))**4)*yo(ii)
	ap(4,8)=ap(4,8)+(((xo(ii))**3)*yo(ii)*yo(ii))
	ap(5,5)=ap(5,5)+(xo(ii)*xo(ii)*yo(ii)*yo(ii))
	ap(5,6)=ap(5,6)+(xo(ii)*yo(ii)*yo(ii)*yo(ii))
	ap(5,7)=ap(5,7)+(((xo(ii))**3)*yo(ii)*yo(ii))
	ap(5,8)=ap(5,8)+(xo(ii)*xo(ii)*((yo(ii))**3))
	ap(6,6)=ap(6,6)+(yo(ii)*yo(ii)*yo(ii)*yo(ii))
	ap(6,7)=ap(6,7)+(xo(ii)*xo(ii)*((yo(ii))**3))
	ap(6,8)=ap(6,8)+(xo(ii)*((yo(ii))**4))
	ap(7,7)=ap(7,7)+(((xo(ii))**4)*yo(ii)*yo(ii))
	ap(7,8)=ap(7,8)+(((xo(ii))**3)*((yo(ii))**3))
	ap(8,8)=ap(8,8)+(((xo(ii))**2)*((yo(ii))**4))
	enddo
c
	do ij=2,8
	do ji=1,ij
	ap(ij,ji)=ap(ji,ij)
	enddo
	enddo
	kk=1
	do ik=1,16
	matp(kk)=1.0d0
	matp(kk+1)=xo(ik)
	matp(kk+2)=yo(ik)
	matp(kk+3)=xo(ik)*xo(ik)
	matp(kk+4)=xo(ik)*yo(ik)
	matp(kk+5)=yo(ik)*yo(ik)
	matp(kk+6)=xo(ik)*xo(ik)*yo(ik)
	matp(kk+7)=xo(ik)*yo(ik)*yo(ik)
	kk=kk+8
	enddo
	call dlinrg(8,ap,8,iap,8)
	is=1
	do ik=1,16
	BMATX(1)=BMATX(1)+matp(is)*mmx(ik)
	BMATX(2)=BMATX(2)+matp(is+1)*mmx(ik)
	BMATX(3)=BMATX(3)+matp(is+2)*mmx(ik)
	BMATX(4)=BMATX(4)+matp(is+3)*mmx(ik)
	BMATX(5)=BMATX(5)+matp(is+4)*mmx(ik)
	BMATX(6)=BMATX(6)+matp(is+5)*mmx(ik)
	BMATX(7)=BMATX(7)+matp(is+6)*mmx(ik)
	BMATX(8)=BMATX(8)+matp(is+7)*mmx(ik)
c
	BMATY(1)=BMATY(1)+matp(is)*mmy(ik)
	BMATY(2)=BMATY(2)+matp(is+1)*mmy(ik)
	BMATY(3)=BMATY(3)+matp(is+2)*mmy(ik)
	BMATY(4)=BMATY(4)+matp(is+3)*mmy(ik)
	BMATY(5)=BMATY(5)+matp(is+4)*mmy(ik)
	BMATY(6)=BMATY(6)+matp(is+5)*mmy(ik)
	BMATY(7)=BMATY(7)+matp(is+6)*mmy(ik)
	BMATY(8)=BMATY(8)+matp(is+7)*mmy(ik)
c
	BMATXY(1)=BMATXY(1)+matp(is)*mmxy(ik)
	BMATXY(2)=BMATXY(2)+matp(is+1)*mmxy(ik)
	BMATXY(3)=BMATXY(3)+matp(is+2)*mmxy(ik)
	BMATXY(4)=BMATXY(4)+matp(is+3)*mmxy(ik)
	BMATXY(5)=BMATXY(5)+matp(is+4)*mmxy(ik)
	BMATXY(6)=BMATXY(6)+matp(is+5)*mmxy(ik)
	BMATXY(7)=BMATXY(7)+matp(is+6)*mmxy(ik)
	BMATXY(8)=BMATXY(8)+matp(is+7)*mmxy(ik)
c
	BMATXZ(1)=BMATXZ(1)+matp(is)*stxz(ik)
	BMATXZ(2)=BMATXZ(2)+matp(is+1)*stxz(ik)
	BMATXZ(3)=BMATXZ(3)+matp(is+2)*stxz(ik)
	BMATXZ(4)=BMATXZ(4)+matp(is+3)*stxz(ik)
	BMATXZ(5)=BMATXZ(5)+matp(is+4)*stxz(ik)
	BMATXZ(6)=BMATXZ(6)+matp(is+5)*stxz(ik)
	BMATXZ(7)=BMATXZ(7)+matp(is+6)*stxz(ik)
	BMATXZ(8)=BMATXZ(8)+matp(is+7)*stxz(ik)
c
	BMATYZ(1)=BMATYZ(1)+matp(is)*styz(ik)
	BMATYZ(2)=BMATYZ(2)+matp(is+1)*styz(ik)
	BMATYZ(3)=BMATYZ(3)+matp(is+2)*styz(ik)
	BMATYZ(4)=BMATYZ(4)+matp(is+3)*styz(ik)
	BMATYZ(5)=BMATYZ(5)+matp(is+4)*styz(ik)
	BMATYZ(6)=BMATYZ(6)+matp(is+5)*styz(ik)
	BMATYZ(7)=BMATYZ(7)+matp(is+6)*styz(ik)
	BMATYZ(8)=BMATYZ(8)+matp(is+7)*styz(ik)
c
	is=is+8
	enddo
c
	do jj=1,8
	do ij=1,8
	MATAX(jj,kj)=MATAX(jj,kj)+iap(jj,ij)*BMATX(ij)
	MATAY(jj,kj)=MATAY(jj,kj)+iap(jj,ij)*BMATY(ij)
	MATAXY(jj,kj)=MATAXY(jj,kj)+iap(jj,ij)*BMATXY(ij)
	MATAXZ(jj,kj)=MATAXZ(jj,kj)+iap(jj,ij)*BMATXZ(ij)
	MATAYZ(jj,kj)=MATAYZ(jj,kj)+iap(jj,ij)*BMATYZ(ij)
	enddo
	enddo
c
	return
	end