c ..................................................................
c
c  mq1dlsc - a multi-quadric radial basis function method
c            with Least Square method
c  Version "c" - the alfa-coefficients and C_j are to be found
c
c           test 1: 1D function interpolation
c           test 2: 1D elliptic equation
c
c ..................................................................
c
	program mq1dlsa
c
	IMPLICIT REAL*8(A-H,O-Z)
c
	parameter (NMAX=501, NXTM=1001, LWA=(NMAX*(3*NMAX+13)/2+1))
	character *10 par
	character *30 comnd
c
c REM	ASA(NMAX,NMAX)
c
	dimension        y(NMAX),          F(NMAX),     C(NMAX),
     *    A(NMAX,NMAX),  FVEC(NMAX),       W(NMAX),  IPVT(NMAX),
     *    WA(LWA)
c
	common/mq1mc/M,NN,ASA(NMAX,NMAX),AAA(NMAX,NMAX),B(NMAX)
	common/mq2pp/ifun,mode,NXT,NQQ,hxt,x(NMAX),xt(NXTM)
c
     	data ifun/1/, mode/1/
c
	external hlsrec
c
c declaration operator-functions
c
	fun(z)   =ffun(ifun,z)
	dfdx(z)  =dfundx(ifun,z)
	d2fdx2(z)=d2fundx(ifun,z)
c
	g(j,z) =sqrt((z-x(j))**2+C(j)**2)

	gx(j,z) =(z-x(j))/sqrt((z-x(j))**2+C(j)**2)

	g2x(j,z)=     1.0/sqrt((z-x(j))**2+C(j)**2) -
     -      (z-x(j))**2/sqrt(((z-x(j))**2+C(j)**2)**3)
c
c  Input solver parameters
c
	in=5
	iout=6
	io=9
	io2=10
	open(io, file='mqle.dat',status='unknown')
cccc	open(io2,file='mqle2.dat',status='unknown')

  1	write(iout,700)
700	format(' mq1dlse:'/
     *	' Choose mode: 1 - 1D interpol. 2 - 1D diff. eq: ',$)
	read(in,*,end=77,err=77)mode
	if((mode.ne.1).and.(mode.ne.2))goto 1
c
  2	write(iout,800)
800	format(' Choose function: 1=x**2, 2=x*(x-1)**2, 3=sin(2*pi*x),'/
     *  '4=exp(-x), 5=abs(x-0.5), 6=atan(10*(x-0.5)): ',$)
     	read(in,*,end=1,err=1)ifun

  3	write(iout,850)
850	format(' Input N, NXT: ',$)
	read(in,*,end=2,err=2)N,NNXT
        NXT=IABS(NNXT)
        NXT=(NXT-1)/2*2+1
        M=N
c

  4	write(iout,900)
900	format(' Choose K (K=1: c(j)=Cmin*(j/N)^(3/4);'/
     *  10x,'(K=2: c(j)=Cmin*((Cmax/Cmin)^((j-1)/(N-1)))'/
     *	' Input K, Cmin and Cmax: ',$)
	read(in,*,end=3,err=3)K,CMIN,CMAX
c
	write(iout,950)N,CMIN,CMAX
950	format(' Used: N =',I5,' Cmin Cmax: ',1p2e15.5)
c
c Generate Cj
c
	do J=1,N
	 if(K.eq.1)then
	  if(CMAX.GT.0.0)then
	    C(J)=CMAX*(FLOAT(J)/FLOAT(N))**0.75
	  else
	    C(J)=ABS(CMAX)
	  endif
	 else
	  C(J)=CMIN*(CMAX/CMIN)**(FLOAT(J-1)/FLOAT(N-1))	  
	 endif
	end do
c
	write(iout,500)' SHAPE PARAMETERS',(I,C(I),I=1,N)
500	format(1x,a/(5(i3,1pe12.4)))
c
c   Set integration points
c
	xt1=0
	xt2=1
	xt(1)=xt1
	dxt=(xt2-xt1)/(NXT-1)
	hxt=dxt
	do I=2,NXT
	 xt(I)=xt(I-1)+dxt
	end do
c
	write(iout,580)' INTEGRATION POINTS',(I,xt(I),I=1,N)
580	format(1x,a/(5(i4,1pe12.3)))
c
c
c   Set Xj points - function knots
c
	x1=0
	x2=1
	x(1)=x1
	dx=(x2-x1)/(N-1)
	do I=2,N
	 x(I)=x(I-1)+dx
	end do
c
	write(iout,600)' FUNCTION POINTS',(I,x(I),I=1,N)
600	format(1x,a/(5(i3,1pe12.4)))
c
c  Select case : mode.eq.1 - interpolation
c                mode.eq.2 - elliptic equation
c
	goto (610,620), mode
c
c  Case  mode.eq.1 - interpolation
c
c  Form a matrix
c
610	continue
c	NN=N+N
	NN=N
	do J=1,NN
	do I=1,NN
	 A(I,J)=0.0
	end do
	end do
c

	do J=1,N
	do I=1,N
	   A(I,J)=g(j,x(i))
	end do
	end do
c
c  Form a RHS
c
	do I=1,N
	   B(I)=FUN(x(i))
	end do
c .....................................................................
c
c   Stiff B.C. (function at the boundary)
c
c .....................................................................
c
      if(NNXT.LT.0)then
c
c   Set B.C.
c
	I=1
	do J=1,N
	 A(I,J)=g(j,x(i))
	 ASA(I,J)=A(I,J)
	end do
	B(I)=FUN(x(i))
	
	I=N
	do J=1,N
	 A(I,J)=g(j,x(i))
	 ASA(I,J)=A(I,J)
	end do
	B(I)=FUN(x(i))
      endif
	goto 650
c
c   End of interpolation - matrix is done.
c
c
c  Case  mode.eq.2 - elliptic equation Uxx=f
c
c  Form a matrix
c	
620	pause 'Mode 2 is not ready yet. Stop.'
c
	do J=1,N
	do I=1,N
	A(I,J)=g2x(j,x(i))
	ASA(I,J)=A(I,J)
	end do
	end do
c
c  Form a RHS
c
	do I=1,N
	 B(I)=d2fdx2(x(i))
	end do
c
c   Set B.C.
c
	I=1
	do J=1,N
	 A(I,J)=g(j,x(i))
	 ASA(I,J)=A(I,J)
	end do
	B(I)=FUN(x(i))
	
	I=N
	do J=1,N
	 A(I,J)=g(j,x(i))
	 ASA(I,J)=A(I,J)
	end do
	B(I)=FUN(x(i))
c
c   End of PDE - matrix is done.
c

c
c Decomposition
c
650	continue
	call decomp(NMAX,N,A,COND,IPVT,W)
	write(iout,1000) N,NXT, COND
1000	format(' MQLS: knots',I5,' points',I5, ' COND =',1pe15.5)
c
c Solve by HMIN/HYBRID
c
	do I=1,N
	Y(I)=B(I)
	end do

	call solve(NMAX,N,A,Y,IPVT)
	write(iout,600)' SOLUTION PRELIMINARY alfa(j)',(I,Y(I),I=1,N)
c
	tol=1.d-12
	info=0
	do I=1,N
	 y(I)=C(I)
	end do
c
	write(iout,*)' Call HYBR1: NN,tol,info,lwa =',NN,tol,info,lwa 
	call  hybrd1(hlsrec,NN,y,fvec,tol,info,wa,lwa)

	write(iout,2000)tol,info
2000	format(' Solved O.K.',' TOL req =',1PE12.3,' INFO =',I5)
c
	goto (201,202,203,204,205),info+1
201	write(iout,*)'Improper input parameters.'
	goto 210

202	write(iout,*)'Relative error between x and the solution ',
     *	             'is at most TOL.'
	goto 210

203	write(iout,*)'Number of calls to fcn has reached or exceeded ',
     *	             '200*(N+1).'
	goto 210

204	write(iout,*)'TOL is too small. no further improvement in ',
     *	             'the approximate solution x is possible.'
	goto 210

205	write(iout,*)'Iteration is not making good progress.'
	goto 210
c
c     OPTMIZATON COMPLETED
c
210	continue
c	write(iout,600)' SOLUTION COEFFICIENTS alfa(j)',(I,Y(I),I=1,N)
	write(iout,600)' DISTANSE L(j)',(I,Y(I),I=1,N)
	do I=1,N
	 C(I)=Y(I)
	end do
c
c   Solve again using new C
c

	do J=1,N
	do I=1,N
	   A(I,J)=g(j,x(i))
	 ASA(I,J)=A(I,J)
	end do
	end do
c
	do I=1,N
	   B(I)=FUN(x(i))
	end do
c
	call decomp(NMAX,N,A,COND,IPVT,W)
	write(iout,1000) N,NXT,COND
	do I=1,N
	Y(I)=B(I)
	end do

	call solve(NMAX,N,A,Y,IPVT)
	write(iout,600)' SOLUTION **FINAL** alfa(j)',(I,Y(I),I=1,N)

c
c	Print solution at the nodes
c
ccc	write (iout,3000)N

	rewind io
	write (io  ,3000)N,COND
3000	format('# SOLUTION ',I5,' POINTS, COND=',1pe12.2/
     *	'#',4x,'I',7x,'X',15x,'NUMER',10x,'ANAL',10x,'ERR',10x,
     *  'Ux (NUMER and ANAL)',10x,'d/dxERR')
c
	NPMX=201
	XFUN=0.0
	DXFUN=0.0
	XXF=0.
	XDX=0.
	do K=1,NPMX
	z=x1+(x2-x1)*FLOAT(K-1)/FLOAT(NPMX-1)
	 SOL=0.0
	 DSDX=0.0
	  do J=1,N
	   SOL=SOL+y(J)*g(J,z)
	   DSDX=DSDX+y(J)*gx(J,z)
	  end do
c
	ANAL=FUN(z)
	ERR=ANAL-SOL
	
	ADFDX=DFDX(Z)
	ERRDF=ADFDX-DSDX

	if(ABS(ERR).GT.ABS(XFUN)) then
	 XFUN=ERR
	 XXF=z
	endif

	if(ABS(ERRDF).GT.ABS(DXFUN)) then
	 DXFUN=ERRDF
	 XDX=z
	endif

c
c	write (iout,3100)K,z,SOL,ANAL,ERR,DSDX,ADFDX,ERRDF
	write (io  ,3100)K,z,SOL,ANAL,ERR,DSDX,ADFDX,ERRDF
3100	format(I5,1P7E15.5)
	end do
	rewind io
c
	write (iout,3150)XFUN,XXF,DXFUN,XDX
3150	format(' Max. Func. error ',1pe15.5,' at ',1pe15.5/
     *	       ' Max. dF/dx error ',1pe15.5,' at ',1pe15.5)
c
c	write solution at the nodes
c
	open(io2,file='mqle2.dat',status='unknown')
	rewind io2
	write (io2  ,3060)N,COND,XFUN,XXF
3060	format('# SOLUTION ',I5,' POINTS, COND=',1pe12.2,
     *         ' Max F,Fx error=',1p2e12.2/
     *	'#',4x,'I',7x,'X',15x,'NUMER',10x,'ANAL',10x,'ERR',10x,
     *  'Ux (NUMER and ANAL)',10x,'d/dxERR'
     *  ,10x,'C[J]',10x,'Curv',10x,'Rad',12x,'Uxx')
c
	XFUN=0.0
	DXFUN=0.0
	XXF=0.
	XDX=0.
	do K=1,N
	z=x(K)
	 SOL=0.0
	 DSDX=0.0
	  do J=1,N
	   SOL=SOL+y(J)*g(J,z)
	   DSDX=DSDX+y(J)*gx(J,z)
	  end do
c
	ANAL=FUN(z)
	ERR=ANAL-SOL
	
	ADFDX=DFDX(Z)
	ERRDF=ADFDX-DSDX

	if(ABS(ERR).GT.ABS(XFUN)) then
	 XFUN=ERR
	 XXF=z
	endif

	if(ABS(ERRDF).GT.ABS(DXFUN)) then
	 DXFUN=ERRDF
	 XDX=z
	endif
	CK=ABS(C(K)/Y(K))
	UXX=d2fdx2(z)
	RCUR=UXX/(SQRT(1.0+ADFDX**2))**3
	RADC=300.
	IF(ABS(RCUR).GT.1./RADC)
     *	RADC=1.0/ABS(RCUR)
c
	write (io2 ,3250)K,z,SOL,ANAL,ERR,DSDX,ADFDX,ERRDF,CK,RCUR,RADC,
     *	                 UXX      
3250	format(I5,1P11E15.5)
	end do
	rewind io2
	close(io2)
c ..........................................................................
c
c   Test points
c
  20	write (iout,3500)
3500	format(' Input X (or q - end, g - gnuplot, s - save): ',$)
	read(in,'(a)',err=66)par
	if(par(1:1).eq."Q")goto 77
	if(par(1:1).eq."q")goto 66
	if(par(1:1).eq."s")then
	 write(iout,3580)
3580	 format(' File name to copy to: ',$)
	 read(in,'(a)',err=66)par
	 comnd='cp -p mqle2.dat '//par
	 write(iout,'(a,a)')' comnd = ',comnd
	 call system(comnd)
	 goto 20
	endif
	if(par(1:1).eq."g")then
cccc	 call system('gnuplot')
	 call gnuplot
	 goto 20
	else
	 read(par,*,err=66)Z
	endif
c
c find numerical solution
c
	 SOL=0.0
	 DSDX=0.0
	 do J=1,N
	  SOL=SOL+y(J)*g(J,z)
	   DSDX=DSDX+y(J)*gx(J,z)
	 end do
c
	ANAL=FUN(z)
	ERR=ANAL-SOL
	ADFDX=DFDX(Z)
	ERRDF=ADFDX-DSDX
c
	write (iout,3800)z,SOL,ANAL,ERR,DSDX,ADFDX
3800	format(7x,'X',9x,'NUMER',8x,'ANAL',9x,'ERR',
     *	10x,'dU/dx (NUMER and ANAL)'/1P6E13.5)

	goto 20
66 	continue
c
	goto 4
c
77	stop
	end
	function ffun(I,x)
c
c  ffun - choice of 5 functions
c
	IMPLICIT REAL*8(A-H,O-Z)
c
	goto (1,2,3,4,5,6),I
1	ffun=x**2
	return
2	ffun=x*(x-1)**2
	return
3	ffun=sin(2*3.14159265358979323844*x)
	return
4	ffun=exp(-x)
	return
5	ffun=abs(x-0.5)
	return
6	ffun=atan(10.0*(x-0.5))
	return
	end
	function dfundx(I,x)
c
c  dfun - derivative: choice of 5 functions
c
	IMPLICIT REAL*8(A-H,O-Z)
c
	goto (1,2,3,4,5,6),I
1	dfundx=2*x
	return
2	dfundx=(x-1)**2 + x*2*(x-1)
	return
3	pi=3.14159265358979323844
	pi2=2.0*pi
	dfundx=pi2*cos(pi2*x)
	return
4	dfundx=-exp(-x)
	return
5	dfundx=-1
	if(x.ge.0.5)dfundx=1
	return
c6	dfundx=10.0/(1.0+tan(10.0*(x-0.5))**2)
c6	dfundx=10.0/(1.0+(x-0.5)**2)
6	del=0.001
	dfun=atan(10.0*(x+del-0.5))-atan(10.0*(x-del-0.5))
	dfundx=dfun/(2.0*del)
	return
	end
	function d2fundx(I,x)
c
c  d2fundx - derivative: choice of 5 functions
c
	IMPLICIT REAL*8(A-H,O-Z)
c
	goto (1,2,3,4,5,6),I
1	d2fundx=2
	return
2	d2fundx=2*(x-1) + 2*(x-1) + x*2
	return
3	pi=3.1415926535897932384
	pi2=2.0*pi
	d2fundx=-pi2**2*sin(pi2*x)
	return
4	d2fundx=exp(-x)
	return
5	d2fundx=0.0
	if(x.eq.0.5)d2fundx=1e6
	return
6	del=0.001
	d2  =atan(10.0*(x+del-0.5))-2*atan(10.0*(x-0.5))
     +	    +atan(10.0*(x-del-0.5))
	d2fundx=d2/del**2
	return
	end
	subroutine gnuplot
c
c   gnuplot - plot the solution, error, derivatives etc
c
	lgp=12
	open(lgp,file='mqle.gp', status='unknown')
	write(lgp,1000)
1000	format('set title "MQ-LSH METHOD INTERPOLATION ERRORS"'/
     * 'set xlabel "X"'/'set ylabel "F, dF/dx ERROR"'/
     * 'plot "mqle.dat" u 2:5 w l 1,"mqle.dat" u 2:8 w l -1'/'pause -1'/
     * 'set title "MQ-LS : F - ORIGINAL AND INTERPOLATED FUNCTION"'/
     * 'set ylabel "F"'/
     * 'plot "mqle.dat" u 2:3 w l 1,"mqle.dat" u 2:4 w l 3,\\'/
     *      '"mqle2.dat" u 2:3 w p 3'/'pause -1'/
     * 'set title "MQ-LS : dF/dx - ORIGINAL AND INTERPOLATED FUNCTION"'/
     * 'set ylabel "dF/dx"'/
     * 'plot "mqle.dat" u 2:6 w l 1,"mqle.dat" u 2:7 w l 5,\\'/
     *      '"mqle2.dat" u 2:6 w p 3'/'pause -1'/
     * 'q')
c
	close (lgp)
c
	write(6,*)'Starting GNUPLOT. Hit RETURN for new figure'
     	 call system('gnuplot mqle.gp')
     	 return
     	 end
	subroutine mxmpy(ndim,n,a,x,y)
c
c   mxmpy - matrix multiply by X:  y = A * x
c
c     input..
c
c        ndim = declared row dimension of the array containing  a.
c
c        n = order of the matrix.
c
c        a = matrix to be triangularized.
c	 x - input vector
c
c     output..
c	 y = A * x
c
      integer ndim,n
      double precision a(ndim,n),x(n),y(n)
c
      do i=1,n
       y(i)=0.0
       do j=1,n
        y(i)=y(i)+A(i,j)*x(j)
       end do
      end do
      return
      end
        subroutine hlsrec(na,x,fvec,iflag)
c         ----------
c         calculate the functions at x and
c         return this vector in fvec.
c         ---------
c
	IMPLICIT REAL*8(A-H,O-Z)
c
	parameter (NMAX=501, NXTM=1001)
        integer na,iflag
        double precision x(na),fvec(na)
c
	dimension W(NMAX),  IPVT(NMAX),    Y(NMAX),    C(NMAX)
	common/mq1mc/M,NN,ASA(NMAX,NMAX),AAA(NMAX,NMAX),B(NMAX)
	common/mq2pp/ifun,mode,NXT,NQQ,hxt,xk(NMAX),xt(NXTM)
c
	data nit/0/
c
c declaration operator-functions
c
	fun(z)   =ffun(ifun,z)
	dfdx(z)  =dfundx(ifun,z)
	d2fdx2(z)=d2fundx(ifun,z)

	g(j,z) =sqrt((z-xk(j))**2+C(j)**2)

	gx(j,z) =(z-xk(j))/sqrt((z-xk(j))**2+C(j)**2)
ccc	gxc(j,z) =-C(j)*(z-xk(j))/sqrt(((z-xk(j))**2+C(j)**2)**3)
	gxc(j,z) =(z-xk(j))/sqrt(((z-xk(j))**2+C(j)**2)**3)

	g2x(j,z)=     1.0/sqrt((z-xk(j))**2+C(j)**2) -
     -      (z-xk(j))**2/sqrt(((z-xk(j))**2+C(j)**2)**3)
c
c ---
	IPREC=1
	IPR=0
c
	N=M
	NN=N
	nit=nit+1
c
c  Form a matrix using the current X=[x,c]
c
	do I=1,N
	   C(I)=X(I)
	end do
c
	do J=1,N
	do I=1,N
	 AAA(I,J)=g(j,xk(i))
	end do
	end do
c
	do I=1,N
	   Y(I)=FUN(xk(i))
	end do
c
c Preconditioning
c
      if(IPREC.ne.0)then
	do I=1,N
	 B(I)=AAA(I,I)
	end do
c
	do J=1,N
	do I=1,N
	 AAA(I,J)=AAA(I,J)/B(I)
	end do
	end do
c	
	do I=1,N
	   Y(I)=Y(I)/B(I)
	end do
      endif
c
	call decomp(NMAX,N,AAA,COND,IPVT,W)
c	
	if(nit.le.20)
     *	write(45,1300) N,NXT,COND
1300	format(' HSOLVC MQLS: knots',I5,' points',I5, ' COND =',1pe15.5)

	call solve(NMAX,N,AAA,Y,IPVT)
c
c  C_j related matrix block
c
	do I=1,N
	 fvec(i)=0.0

	  do k=1,N
	  uxk=0.0
	  do J=1,N
	   uxk=uxk+y(J)*gx(j,xk(k))
	  end do
c
	   fvec(i)=fvec(i)+Y(I)*gxc(i,xk(k))
     *	           *(uxk-dfdx(xk(k)))
	  end do
c
	end do
c
	resid=0.0
	xx2=0.0

	do i=1,n
	 resid=resid+fvec(i)**2
	 xx2=xx2+x(i)*x(i)
	end do
	resid=sqrt(resid)
	xx2=sqrt(xx2)
c
	if((nit.le.10).or.(nit/10*10.eq.nit))
     *	write(6,100)nit,na,M,xx2,resid,COND
100   format('hlsrec: nit,nn,M=',3I4,' |X| =',1pe11.3,' RESID=',1pe10.2
     *	,' COND=',1pe10.2)
c
        return
        end