! sp1 Geoff's brain production
! Copyright (C) 1999,2010 Gauthier Delerce
!
! This program is free software; you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation; either version 2, or (at your option)
! any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with this program; if not, write to the Free Software Foundation,
! Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
!
! Gauthier Delerce <gauthier@delerce.fr>
subroutine sp1(a,x,y,n,ro,s,t,r,b,aa,c1,c2,c3,max)
implicit real(8) (a-g,o-z)
real x(n),y(n),a(n)
real(8) s(n),t(3),ro
real(8) r(max*(max-1)/2),b(max),aa((max-3)*(max-2)/2)
real(8) c1(max),c2(max),c3(max)
!*----------- computation of kernel matrix K (I,J) I<J -------------------
! write(*,*) ' DEBUT SP1'
! write(*,*) 'max, aa size,n3',max,(max-3)*(max-2)/2,n-3
!c pause
! print*,111,a
! print*,''
! print*,x;print*,''
! print*,y;print*,''
! print*,n;print*,''
! print*,ro;print*,''
! print*,s;print*,''
! print*,t;print*,''
! print*,r;print*,''
! print*,b;print*,''
! print*,aa;print*,''
! print*,c1;print*,''
! print*,c2;print*,''
! print*,c3;print*,''
! print*,max;print*,''
n3=n-3
do j=2,n
jm1=j-1
ij=(jm1*(j-2))/2
do i=1,jm1
dx=dble(x(i)-x(j))
dy=dble(y(i)-y(j))
a3=dx*dx+dy*dy
r(ij+i)=a3*dlog(a3)
!c write(*,*)i,j,ij+i
!c pause
end do
end do
!c pause
!*------------------- computation of C1,C2,C3(I) --------------------------
v1=dble(y(3)-y(1))
v4=dble(x(2)-x(1))
v2=dble(x(1)-x(3))
v3=dble(y(1)-y(2))
a1=v1*v4-v2*v3
v1=v1/a1
v2=v2/a1
v3=v3/a1
v4=v4/a1
do i=1,n3
ip3=i+3
a1=dble(x(1)-x(ip3))
a2=dble(y(1)-y(ip3))
c3(i)=v3*a1+v4*a2
c2(i)=v1*a1+v2*a2
c1(i)=-(1.d0+c2(i)+c3(i))
end do
!c pause
!*-------------------- computation of AA(I,J) I <= J ----------------------
do i=1,n3
i1=(i*(i+1))/2
ij=i1+i+2
a1=c2(i)*r(1)+c3(i)*r(2)+r(ij)
a2=c3(i)*r(3)+r(ij+1)
a3=ro*(c1(i)*c1(i)+c2(i)*c2(i)+c3(i)*c3(i)+1.d0)
aa(i1)=2.d0*(c1(i)*a1+c2(i)*a2+c3(i)*r(ij+2))+a3
end do
do j=2,n3
jm1=j-1
ij=(j*(j+3)+4)/2
it=(j*jm1)/2
do i=1,jm1
i1=it+i
j1=(i*(i+3)+4)/2
a1=c2(j)*r(1)+c3(j)*r(2)+r(ij)
a2=c1(j)*r(1)+c3(j)*r(3)+r(ij+1)
a3=c1(j)*r(2)+c2(j)*r(3)+r(ij+2)
a4=c1(j)*r(j1)+c2(j)*r(j1+1)+c3(j)*r(j1+2)+r(ij+i+2)
a5=ro*(c1(i)*c1(j)+c2(i)*c2(j)+c3(i)*c3(j))
aa(i1)=c1(i)*a1+c2(i)*a2+c3(i)*a3+a4+a5
!c write(*,*) aa(i1)
end do
end do
!c pause
!*------------------ Cholesky decomposition -------------------------------
! print*,1111
call ludecp(aa,aa,n3,d1,d2,ier,(max-3)*(max-2)/2)
!c pause
!*---------- forward/back substitution for A ------------------------------
! print*,2222
do i=1,n3
b(i)=c1(i)*a(1)+c2(i)*a(2)+c3(i)*a(3)+a(i+3)
end do
call luelmp(aa,b,n3,b,(max-3)*(max-2)/2)
do i=4,n
s(i)=b(i-3)
end do
! print*,3333
!*-------------------- compute A, spline coefficients ---------------------
a1=0.d0
a2=0.d0
a3=0.d0
do i=1,n3
ip3=i+3
a1=a1+c1(i)*s(ip3)
a2=a2+c2(i)*s(ip3)
a3=a3+c3(i)*s(ip3)
end do
s(1)=a1
s(2)=a2
s(3)=a3
a1=s(2)*r(1)+s(3)*r(2)+ro*s(1)
a2=s(1)*r(1)+s(3)*r(3)+ro*s(2)
a3=s(1)*r(2)+s(2)*r(3)+ro*s(3)
do i=4,n
i1=(i*(i-3)+4)/2
a1=a1+s(i)*r(i1)
a2=a2+s(i)*r(i1+1)
a3=a3+s(i)*r(i1+2)
end do
a1=a(1)-a1
a2=a(2)-(a2+a1)
a3=a(3)-(a3+a1)
t(1)=v1*a2+v3*a3
t(2)=v2*a2+v4*a3
t(3)=a1-(x(1)*t(1)+y(1)*t(2))
! print*,7777
!* ----------- return with spline coefficients in S & T ------------------
return
end