!MS$real:8	
	program START

	use AVDef    
	use DFLIB
	
	integer, parameter:: dil=150,dl=64,del=31 !150

      integer, parameter::n=(dil)*4,k=n/2,nl=n/2,np=n/2
	  	common Ph0,Ph1,Ph2,Ph3,Ph4,Pha0T,oleole
	Common PhN1,PhN2,PhN3,PhN1T
	real z(N*(N-K+1)+(N-K+1)*(N-K+2)/2+1)
     real r1(N*(N-K+1)),z1(N*(N-K+1)),y(n+1),a1(N*(N+1)),oleole(N*(N+1))
     real om((N-K+1)*(N-K+2)/2)
     real rew(dil+4,31),ww(dl,dil+4),ww1(dil+4,dl)
	 real Aleft(K*(n+1)),Aright(K*(N+1)),AMat1(n/2,n+1),Amat2(n/2,n+1)
	 dimension itp(n)
	integer num1
	real xx(-5:dil+4)
	real ess(dil),x,ds,esd,ess1(n/4),ess2(n/4)
	real Mm(del,dl)
      real sina(dl)
		real suk,vnu,tempv
		character key
	real Pha0T(n/4,n/4)
	real Ph0(n/4,n/4),Ph1(n/4,n/4)
	real Ph2(n/4,n/4),Ph3(n/4,n/4),Ph4(n/4,n/4)
	real PhN1(n/4,n/4),PhN1T(n/4,n/4)
	real PhN2(n/4,n/4),PhN3(n/4,n/4),PhNN(n/4,n/4)
	open (6,file='dior.txt',status='unknown')
	open (7,file='dior1.txt',status='unknown')
		open (8,file='dior2.txt',status='unknown')

		

 itp=(/(i,i=1,n)/)	
!======        ========
! do i=n/3 + 1,n/2
!	itp(i) = itp(i) + n/6;
! enddo

! do i=n/2 + 1,2*n/3
!	itp(i) = itp(i) - n/6;
! enddo
!======================================================================

!======    =======================================
!do i=n/6 + 1,n/3
!  itp(i) = itp(i) + 4*n/6;
!enddo
 
!do i=5*n/6 + 1,n
!  itp(i) = itp(i) - 4*n/6;
!enddo

!AMat1=0.; AMat2=0.;

! do i=1,n/6
! AMat1(i,i)=1.
! enddo

! do i=n/6 + 1,n/3
! AMat1(i+n/6,i)=1.
! enddo
  
! do i=n/3 + 1,n/2
! AMat1(i-n/6,i)=1.
! enddo

!------------------

! do i=1,n/6
! AMat2(i,i)=1.
! enddo

! do i=5*n/6 + 1,n
! AMat2(i-n/2,i)=1.
! enddo
  
! do i=n/3 + 1,n/2
! AMat2(i-n/6,i)=1.
! enddo





!======================================================================

  

!Aleft=0.
!Aright=0.

!	  x=0.
!	esd=10.

!	xx(0)=0.
!	ds=esd/(dil-1)
!	do i=1,dil+4
!	xx(i)=xx(i-1)+ds
!	enddo
!	i=-1
!	do while (i>-6)
!	xx(i)=xx(i+1)-ds
!	i=i-1
!	enddo
!	do i=0,dil-2
!	if (mod(i,2)==0.and.i/=1) then
!	ess(i+1)=xx(i)+(1./2.-sqrt(3.)/6.)*ds
!	ess(i+2)=xx(i)+(1./2.+sqrt(3.)/6.)*ds
!	end if 
!	enddo

!	do j=1,n/4
!	do i=1,n/4
!	Ph0(i,j)=phi(0,j-1,dil-1,xx,ess(i))
!	Ph1(i,j)=phi(1,j-1,dil-1,xx,ess(i))
!	Ph2(i,j)=phi(2,j-1,dil-1,xx,ess(i))
!	Ph3(i,j)=phi(3,j-1,dil-1,xx,ess(i))
!	Ph4(i,j)=phi(4,j-1,dil-1,xx,ess(i))
!	enddo
!	enddo

!vnu=.3
!do j=1,n+1; do i=1,n/2
!Aleft((j-1)*N/2+i)=AMat1(i,j)
!Aright((j-1)*N/2+i)=AMat2(i,j)
! enddo
! enddo


!===================================================
AMat1=0.; AMat2=0.;

Aleft=0.
Aright=0.
!--     ------------------------------
open(102,file='../BqLeft.txt',status='unknown')
      do i=1,n/2
	  do j=1,n 
		read(102,*)AMat1(i,j)
		
	  enddo
   enddo
close(102)

open(102,file='../Bq.txt',status='unknown')
      do i=1,n/2
	  do j=1,n 
		read(102,*)AMat2(i,j)
		
	  enddo
   enddo
close(102)

!do j=1,n 
!	do i=1,n/2
!		if (i==j) then
!			Aright((j-1)*N/2+i) = 1
!			Aleft((j-1)*N/2+i) = 1
!		end if
!	enddo
!enddo


!--      . Left ~ z = 0,  Right ~ z = c.
open(103,file='../f0.txt',status='unknown')
	  do i=1,n/2
		read(103,*)AMat1(i,n + 1)
	  enddo
close(103)

open(104,file='../fc.txt',status='unknown')
	  do i=1,n/2
		read(104,*)AMat2(i,n + 1)
	  enddo
close(104)

do j=1,n+1 
do i=1,n/2
Aleft((j-1)*N/2+i)=AMat1(i,j)
Aright((j-1)*N/2+i)=AMat2(i,j)
enddo
enddo
!===================================================

!--   --------------------------------------
	open(105,file='../A.txt',status='unknown')
      do j=1,n 
	  do i=1,n
		read(105,*)oleole((j-1)*N+i)
	  enddo
	  enddo
	close(105)

    !--     -------------------------
	open(106,file='../f.txt',status='unknown')
	  do i=1,n
		read(106,*)oleole(N*N+i)
	  enddo
	close(106)
!---------------------------------------------------------
!=========================================================
!---------------------------
open (666,file='Aleft.txt',status='unknown')
write (666,1001)Aleft
close(666)

open (667,file='Aright.txt',status='unknown')
write (667,1001)Aright
close(667)

1001  FORMAT(1(E12.5,2X))
!---------------------------

PhN1=0.;PhN2=0.;PhN3=0.

!	call AlterOb(Ph0,Pha0T,n/4)	

	suk=esd/(dl-1)
	do i=0,dl-1
		sina(i+1)=i*suk
	enddo

      N1=(N-K+1)*N
      M1=(N-K+1)*(N-K+2)/2
      NZ=N1+M1+1
      NY=N+1
      NN=N*NY
      CALL DO(N,K,X,H,N1,M1,NZ,NY,NN,NL,NP, &
     Z,R1,Z1,Y,A1,OM,Aleft,Aright,ITP)
  
  close(6)

!  open(6,file='dior.txt',status='unknown')
	
!	kk=del
!	do while (kk>0)
!read(6,*)ww(kk,1:n/4)
!	kk=kk-1
!	enddo

!	do l=1,dl
!	do kk=1,del
!mm(kk,l)=0.
!do i=1,n/4
!mm(kk,l)=mm(kk,l)+ww(kk,i)*phi(0,i-1,dil-1,xx,sina(l))
!	enddo
! 	enddo;		enddo

!	 call faglStartWatch(Mm, status)

!	print *, "Starting Array Viewer"

!	call faglShow(Mm, status)
!		key = GETCHARQQ()
!	call faglEndWatch(Mm, status)        

print *, "Complete"
key = GETCHARQQ()
 
END


	  SUBROUTINE FORWAR(N,M,X,H,IP1,IP2,IP3,NA,A)
      real(8) A(NA),h,x
	  IP1=1; IP2=0;   IP3=1
      
	  H=0.01

	  IF(ABS(0.-X).LE.H/2.) IP2=1
	  IF(ABS(0.1-X).LE.H/2.) IP2=1
	  IF(ABS(0.2-X).LE.H/2.) IP2=1
	  IF(ABS(0.3-X).LE.H/2.) IP2=1
	  IF(ABS(0.4-X).LE.H/2.) IP2=1
	  IF(ABS(0.5-X).LE.H/2.) IP2=1
	  IF(ABS(0.6-X).LE.H/2.) IP2=1
	  IF(ABS(0.7-X).LE.H/2.) IP2=1
	  IF(ABS(0.8-X).LE.H/2.) IP2=1
	  IF(ABS(0.9-X).LE.H/2.) IP2=1
	  IF(ABS(1.-X).LE.H/2.) IP2=1

      IF(ABS(1.-X).LE.H/2.) IP3=0
      RETURN
      END

      SUBROUTINE RESULT(N,NI,A)
      real(8) A(N)
		!write (6,1)a(2:n)
		do i=2,n
		   write (6,1)a(i)
		1  FORMAT(600(E12.5,1X))
		enddo

      RETURN
      END
 
!HERE MUST DEFINE: Pha0T, Ph4, Pha0T, Ph2, qf - right side vector
     
    SUBROUTINE MATRIX(NN,X,NA,A)
	!MS$real:8
	integer, parameter:: dil=54,n=(dil)*4
	common Ph0,Ph1,Ph2,Ph3,Ph4,Pha0T,oleole
	Common PhN1,PhN2,PhN3,PhN1T
      real A(NA),oleole(N*(N+1))
	real h,x,vn_st,vn_ts !coef  of  Puasson
	real tmp4(n/4,n/4)
	real qdz(n/4),qqdz(n/4),qf(n/4)
	real dd(4,4,n/4,n/4),dd1((n/4)*4,(n/4)*4)
	real Pha0T(n/4,n/4),Phb0T(n/4,n/4)
	real Ph0(n/4,n/4),Ph1(n/4,n/4)
	real Ph2(n/4,n/4),Ph3(n/4,n/4),Ph4(n/4,n/4)
	real PhN1(n/4,n/4),PhN1T(n/4,n/4)
	real PhN2(n/4,n/4),PhN3(n/4,n/4)
	real Test(n/4 * n/4)
!===========================================
	real E1, E2, G12, niu1, niu2, D11, D12, D22, D66 , hhh
	real aaa, bbb, qqq

!	hhh = 0.01
!	E1 = 4.76 * 10000
!	E2 = 2.07 * 10000
!	G12 = 0.531 * 10000
!	niu1 = 0.149
!	niu2 = 0.0647

!	D11 = hhh * E1/(1-niu1*niu2)
!	D12 = hhh * niu1 * E2/(1-niu1*niu2)
!	D22 = hhh * E2/(1-niu1*niu2)
!	D66 = hhh * G12

!	aaa = -D22/D11
!	bbb = -2*(D12 + 2*D66)/D11

!	qqq = 0.01

!===========================================

!	dd=0.;dd1=0.
!	a=0.
!	do i=1,n/4
!	dd(1,2,i,i)=1.
!	dd(3,4,i,i)=1.
!	dd(2,3,i,i)=1.
!	enddo

!	dd(4,1,:,:)=aaa * matmul(Pha0T,Ph4) !-matmul(Pha0T,Ph4)
!	dd(4,3,:,:)=bbb * matmul(Pha0T,Ph2) !-2.*matmul(Pha0T,Ph2)


!	do j=0,3;	do l=0,3
!	do m=1,n/4
!	do i=1,n/4
!	dd1(j*(n/4)+m,l*(n/4)+i)=dd(j+1,l+1,m,i)
!	enddo;	enddo
!	enddo;	enddo

!	do j=1,(n/4)*4;	do i=1,(n/4)*4
!	a((j-1)*(n/4)*4+i)=dd1(j,i)
!	enddo;				enddo


!-----Right side -----------------------------------------
	
!	do i=1,n/4
!	qdz(i)=qqq/D11 
!	enddo
	 
!	call multmv(n/4,Pha0T,qdz,qf)
	
!	do i=(n/4)*3+1,n
!	a(n*n+i)=qf(i-(n/4)*3)
!	enddo


!--   --------------------------------------
!	open(105,file='../A.txt',status='unknown')
!      do j=1,n 
!	  do i=1,n
!		read(105,*)a((j-1)*N+i)
!	  enddo
!	  enddo
!	close(105)

    !--     -------------------------
!	open(106,file='../f.txt',status='unknown')
!	  do i=1,n
!		read(106,*)a(N*N+i)
!	  enddo
!	close(106)
    !---------------------------------------------------------
a=oleole

!open(106,file='../f.txt',status='unknown')
! do i=1,n
!	read(106,*)a(N*N+i)
! enddo
!close(106)

!do i=1,(n/3)*2
!	a(N*N+i) = a(N*N+i) * Sin(3.14159265*x)
!enddo

!do i=(n/3)*2+1,n
!	a(N*N+i) = a(N*N+i) * Cos(3.14159265*x)
!enddo

!-----------------------------------------------
open (668,file='A.txt',status='unknown')
write (668,1002)A
close(668)

1002  FORMAT(1(E12.5,2X))
!-----------------------------------------------
	
!---------------------------------------------------------
 111  FORMAT(1(E12.5,2X)) 
      RETURN
	 
	  contains
	  subroutine multmv(nu,au,b,c)
	real(8) au(nu,nu),c(nu),b(nu)
	integer i,j,k
	
	do i=1,nu
		c(i)=0.
		do ku=1,nu
			c(i)=c(i)+au(i,ku)*b(ku)
		enddo
	end do
	end subroutine multmv

      END subroutine

!----------------------------------------------------------------------
!--------hold tightly borders -----------------------------------------    
	real(8)	function phi(j,i,n,xx,x)
		real(8) x,xx(-5:n+5)
		integer j,i,n,k
	if (j==0) then	
		if(i==0)then 
	phi=165./4.*splB5(-2,n,xx,x)-33./8.*splB5(-1,n,xx,x)+splB5(0,n,xx,x)
		else if(i==1) then 
	phi=splB5(-1,n,xx,x)-26./33.*splB5(0,n,xx,x)+splB5(1,n,xx,x)
			else if(i==2) then
	phi=splB5(-2,n,xx,x)-1./33.*splB5(0,n,xx,x)+splB5(2,n,xx,x)
				else if(i==n-2) then
	phi=splB5(n+2,n,xx,x)-1./33.*splB5(n,n,xx,x)+splB5(n-2,n,xx,x)
					else if(i==n-1) then
	phi=splB5(n+1,n,xx,x)-26./33.*splB5(n,n,xx,x)+splB5(n-1,n,xx,x)
						elseif (i==n) then	
	phi=165./4.*splB5(n+2,n,xx,x)-33./8.*splB5(n+1,n,xx,x)+splB5(n,n,xx,x)
							else
			phi=splB5(i,n,xx,x)
		endif
	elseif  (j==1)then
	
		if(i==0)then 
	phi=165./4.*splB51(-2,n,xx,x)-33./8.*splB51(-1,n,xx,x)+splB51(0,n,xx,x)
		else if(i==1) then 
	phi=splB51(-1,n,xx,x)-26./33.*splB51(0,n,xx,x)+splB51(1,n,xx,x)
			else if(i==2) then
	phi=splB51(-2,n,xx,x)-1./33.*splB51(0,n,xx,x)+splB51(2,n,xx,x)
				else if(i==n-2) then
	phi=splB51(n+2,n,xx,x)-1./33.*splB51(n,n,xx,x)+splB51(n-2,n,xx,x)
					else if(i==n-1) then
	phi=splB51(n+1,n,xx,x)-26./33.*splB51(n,n,xx,x)+splB51(n-1,n,xx,x)
						elseif (i==n) then	
	phi=165./4.*splB51(n+2,n,xx,x)-33./8.*splB51(n+1,n,xx,x)+splB51(n,n,xx,x)
							else
			phi=splB51(i,n,xx,x)
		endif
	elseif (j==2)then
	
		if(i==0)then 
	phi=165./4.*splB52(-2,n,xx,x)-33./8.*splB52(-1,n,xx,x)+splB52(0,n,xx,x)
		else if(i==1) then 
	phi=splB52(-1,n,xx,x)-26./33.*splB52(0,n,xx,x)+splB52(1,n,xx,x)
			else if(i==2) then
	phi=splB52(-2,n,xx,x)-1./33.*splB52(0,n,xx,x)+splB52(2,n,xx,x)
				else if(i==n-2) then
	phi=splB52(n+2,n,xx,x)-1./33.*splB52(n,n,xx,x)+splB52(n-2,n,xx,x)
					else if(i==n-1) then
	phi=splB52(n+1,n,xx,x)-26./33.*splB52(n,n,xx,x)+splB52(n-1,n,xx,x)
						elseif (i==n) then	
	phi=165./4.*splB52(n+2,n,xx,x)-33./8.*splB52(n+1,n,xx,x)+splB52(n,n,xx,x)
							else
			phi=splB52(i,n,xx,x)
			endif
		elseif (j==3)then
			if(i==0)then 
	phi=165./4.*splB53(-2,n,xx,x)-33./8.*splB53(-1,n,xx,x)+splB53(0,n,xx,x)
		else if(i==1) then 
	phi=splB53(-1,n,xx,x)-26./33.*splB53(0,n,xx,x)+splB53(1,n,xx,x)
			else if(i==2) then
	phi=splB53(-2,n,xx,x)-1./33.*splB53(0,n,xx,x)+splB53(2,n,xx,x)
				else if(i==n-2) then
	phi=splB53(n+2,n,xx,x)-1./33.*splB53(n,n,xx,x)+splB53(n-2,n,xx,x)
					else if(i==n-1) then
	phi=splB53(n+1,n,xx,x)-26./33.*splB53(n,n,xx,x)+splB53(n-1,n,xx,x)
						elseif (i==n) then	
	phi=165./4.*splB53(n+2,n,xx,x)-33./8.*splB53(n+1,n,xx,x)+splB53(n,n,xx,x)
							else
			phi=splB53(i,n,xx,x)
		endif
		elseif (j==4)then
			if(i==0)then 
	phi=165./4.*splB54(-2,n,xx,x)-33./8.*splB54(-1,n,xx,x)+splB54(0,n,xx,x)
		else if(i==1) then 
	phi=splB54(-1,n,xx,x)-26./33.*splB54(0,n,xx,x)+splB54(1,n,xx,x)
			else if(i==2) then
	phi=splB54(-2,n,xx,x)-1./33.*splB54(0,n,xx,x)+splB54(2,n,xx,x)
				else if(i==n-2) then
	phi=splB54(n+2,n,xx,x)-1./33.*splB54(n,n,xx,x)+splB54(n-2,n,xx,x)
					else if(i==n-1) then
	phi=splB54(n+1,n,xx,x)-26./33.*splB54(n,n,xx,x)+splB54(n-1,n,xx,x)
						elseif (i==n) then	
	phi=165./4.*splB54(n+2,n,xx,x)-33./8.*splB54(n+1,n,xx,x)+splB54(n,n,xx,x)
							else
			phi=splB54(i,n,xx,x)
		endif
		else 
		stop
	endif
	return
	end
!---------------------------------------------------------------------- 
!--------on-hinges borders---------------------------------------------
	real(8)	function phi2(j,i,n,xx,x)
		real(8) x,xx(-5:n+5)
		integer j,i,n,k
	if (j==0) then	
		if(i==0)then 
	phi2=165./4.*splB5(-2,n,xx,x)-33./8.*splB5(-1,n,xx,x)+splB5(0,n,xx,x)
		else if(i==1) then 
	phi2=splB5(-1,n,xx,x)-26./33.*splB5(0,n,xx,x)+splB5(1,n,xx,x)
			else if(i==2) then
	phi2=splB5(-2,n,xx,x)-1./33.*splB5(0,n,xx,x)+splB5(2,n,xx,x)
				else if(i==n-2) then
	phi2=-1.*splB5(n+2,n,xx,x)+splB5(n-2,n,xx,x)
					else if(i==n-1) then
	phi2=-1.*splB5(n+1,n,xx,x)+splB5(n-1,n,xx,x)
						elseif (i==n) then	
	phi2=12.*splB5(n+2,n,xx,x)-3.*splB5(n+1,n,xx,x)+splB5(n,n,xx,x)
							else
			phi2=splB5(i,n,xx,x)
		endif
	elseif  (j==1)then
	
		if(i==0)then 
	phi2=165./4.*splB51(-2,n,xx,x)-33./8.*splB51(-1,n,xx,x)+splB51(0,n,xx,x)
		else if(i==1) then 
	phi2=splB51(-1,n,xx,x)-26./33.*splB51(0,n,xx,x)+splB51(1,n,xx,x)
			else if(i==2) then
	phi2=splB51(-2,n,xx,x)-1./33.*splB51(0,n,xx,x)+splB51(2,n,xx,x)
				else if(i==n-2) then
	phi2=-1.*splB51(n+2,n,xx,x)+splB51(n-2,n,xx,x)
					else if(i==n-1) then
	phi2=-1.*splB51(n+1,n,xx,x)+splB51(n-1,n,xx,x)
						elseif (i==n) then	
	phi2=12.*splB51(n+2,n,xx,x)-3.*splB51(n+1,n,xx,x)+splB51(n,n,xx,x)
							else
			phi2=splB51(i,n,xx,x)
		endif
	elseif (j==2)then
	
		if(i==0)then 
	phi2=165./4.*splB52(-2,n,xx,x)-33./8.*splB52(-1,n,xx,x)+splB52(0,n,xx,x)
		else if(i==1) then 
	phi2=splB52(-1,n,xx,x)-26./33.*splB52(0,n,xx,x)+splB52(1,n,xx,x)
			else if(i==2) then
	phi2=splB52(-2,n,xx,x)-1./33.*splB52(0,n,xx,x)+splB52(2,n,xx,x)
				else if(i==n-2) then
	phi2=-1.*splB52(n+2,n,xx,x)+splB52(n-2,n,xx,x)
					else if(i==n-1) then
	phi2=-1.*splB52(n+1,n,xx,x)+splB52(n-1,n,xx,x)
						elseif (i==n) then	
	phi2=12.*splB52(n+2,n,xx,x)-3.*splB52(n+1,n,xx,x)+splB52(n,n,xx,x)
							else
			phi2=splB52(i,n,xx,x)
			endif
		elseif (j==3)then
			if(i==0)then 
	phi2=165./4.*splB53(-2,n,xx,x)-33./8.*splB53(-1,n,xx,x)+splB53(0,n,xx,x)
		else if(i==1) then 
	phi2=splB53(-1,n,xx,x)-26./33.*splB53(0,n,xx,x)+splB53(1,n,xx,x)
			else if(i==2) then
	phi2=splB53(-2,n,xx,x)-1./33.*splB53(0,n,xx,x)+splB53(2,n,xx,x)
				else if(i==n-2) then
	phi2=-1.*splB53(n+2,n,xx,x)+splB53(n-2,n,xx,x)
					else if(i==n-1) then
	phi2=-1.*splB53(n+1,n,xx,x)+splB53(n-1,n,xx,x)
						elseif (i==n) then	
	phi2=12.*splB53(n+2,n,xx,x)-3.*splB53(n+1,n,xx,x)+splB53(n,n,xx,x)
							else
			phi2=splB53(i,n,xx,x)
		endif
		elseif (j==4)then
			if(i==0)then 
	phi2=165./4.*splB54(-2,n,xx,x)-33./8.*splB54(-1,n,xx,x)+splB54(0,n,xx,x)
		else if(i==1) then 
	phi2=splB54(-1,n,xx,x)-26./33.*splB54(0,n,xx,x)+splB54(1,n,xx,x)
			else if(i==2) then
	phi2=splB54(-2,n,xx,x)-1./33.*splB54(0,n,xx,x)+splB54(2,n,xx,x)
				else if(i==n-2) then
	phi2=-1.*splB54(n+2,n,xx,x)+splB54(n-2,n,xx,x)
					else if(i==n-1) then
	phi2=-1.*splB54(n+1,n,xx,x)+splB54(n-1,n,xx,x)
						elseif (i==n) then	
	phi2=12.*splB54(n+2,n,xx,x)-3.*splB54(n+1,n,xx,x)+splB54(n,n,xx,x)
							else
			phi2=splB54(i,n,xx,x)
		endif
		else 
		stop
	endif
	return
	end
