Untitled

program main
 use tt_lib
 use ttaux_lib
 use dmrgg_lib
 use time_lib
 use quad_lib
 use default_lib
 use mat_lib
 use omp_lib
 implicit none
 include 'mpif.h'
 type(dtt) :: tt,qq
 integer :: i,j,p,m,n,nx,r,piv,decay,info,nproc,me,adj
 integer(kind=8) :: neval
 double precision :: f,bnd,t1,t2,tcrs, einf,efro,ainf,afro, acc,val,tru,h,w,t
 double precision,allocatable::par(:),mat(:,:)
 double precision,parameter :: tpi=6.2831853071795864769252867665590057683943387987502116419498891846156328125724179972560696506842341359642961730265646132941876892191011644634507188162569622349005682054038770422111192892458979098607639288576219513318668922569512964675735663305424038182912971338469206972209086532964267872145204982825474491740132126311763497630418419256585081834307287357851807200226610610976409330427682939038830232188661145407315191839061843722347638652235862102370961489247599254991347037715054497824558763660238982596673467248813132861720427898927904494743814043597218874055410784343525863535047693496369353388102640011362542905271216555715426855155792183472743574429368818024499068602930991707421015845593785178470840399122242580439217280688363196272595495426199210374144226999999967459560999021194634656321926371900489189106938166052850446165066893700705238623763420200062756775057731750664167628412343553382946071965069808575109374623191257277647075751875039155637155610643424536132260038557532223918184328403d0
 double precision,parameter :: log2=0.6931471805599453094172321214581765680755001343602552541206800094933936219696947156058633269964186875420014810205706857336855202357581305570326707516350759619307275708283714351903070386238916734711233501153644979552391204751726815749320651555247341395258829504530070953263666426541042391578149520437404303855008019441706416715186447128399681717845469570262716310645461502572074024816377733896385506952606683411372738737229289564935470257626520988596932019650585547647033067936544325476327449512504060694381471046899465062201677204245245296126879465461931651746813926725041038025462596568691441928716082938031727143677826548775664850856740776484514644399404614226031930967354025744460703080960850474866385231381816767514386674766478908814371419854942315199735488037516586127535291661000710535582498794147295092931138971559982056543928717000721808576102523688921324497138932037843935308877482597017155910708823683627589842589185353024363421436706118923678919237231467232172053401649256872747782344535d0
 double precision,parameter :: zeta3=1.202056903159594285399738161511449990764986292340498881792271555341838205786313090186455873609335258146199157795260719418491995998673283213776396837207900161453941782949360066719191575522242494243961563909664103291159095780965514651279918405105715255988015437109781102039827532566787603522336984941661811057014715778639499737523785277937030956025701853182790003076547107563048843320869711573742380793445031607625317714535444411831178182249718526357091824489987962035083357561720226033937858703281312678079900541773486911525370656237057440966221712902627320732361492242913040528555372341033077577798064242024304882815210009146026538220696271552020822743350010152948011986901176259516763669981718355752348807037195557423472940835952088616662025728537558130792825864872821737055661968989526620187768106292008177923381358768284264124324314802821736745067206935076268953043459393750329663637757506247332399234828831077339052768020075798435679371150509005027366047114008533503436467224856531518117766181092d0
 character(len=1) :: a
 character(len=32) :: aa
 logical :: rescale
 double precision,external :: dfunc_ising_discr
 double precision,external :: fun





 tt%l=1;tt%m=6-1;tt%n=27;tt%r=1;call alloc(tt)
  

call dtt_dmrgg(tt,fun,par,maxrank=500,accuracy=1E-7,pivoting=3)

 call dealloc(tt)
 call mpi_finalize(info)
 if(info.ne.0)then;write(*,*)'mpi: finalize fail: ',info;stop;endif
end program

double precision function dfunc_ising_discr(m,ind,n,par) result(f)
 implicit none
 integer,intent(in) :: m
 integer,intent(in) :: ind(m),n(m)
 double precision,intent(inout),optional :: par(*)
 integer :: nodes,weights,i,j,id
 double precision :: uij,a,b,v,w,wk,vk
 id=int(par(2*n(1)+1))
 nodes=0       ! t2 t3 ... tm
 weights=n(1)  ! w2 w3 ... wm
 if(id.eq.2 .or. id.eq.3)then
  a=1.d0
  do i=0,m
   uij=1.d0
   do j=i+1,m
    uij=uij*par(nodes+ind(j))
    a = a * (( uij-1.d0 ) / ( uij+1.d0 ))**2
   end do
  end do
 end if
 if(id.eq.1 .or. id.eq.2)then
  v=1.d0; w=1.d0; vk=1.d0; wk=1.d0
  do i=1,m
   vk=vk*par(nodes+ind(m-i+1))
   wk=wk*par(nodes+ind(i))
   v=v+vk
   w=w+wk
  end do
  b=1.d0/(v*w)
 end if 
 select case(id)
 case(1);f=2*b
 case(2);f=2*a*b
 case(3);f=2*a
 case default
  write(*,*)'unknown id: ',id;stop
 end select 
  
 ! apply weights
 do i=1,m
  f=f*par(weights+ind(i))
 end do
end function




double precision function fun(m, ind, n, par) result(result)
  implicit none
  integer, intent(in) :: m
  integer, intent(in) :: ind(m), n(m)
  double precision, intent(inout), optional :: par(*)
  double precision :: b(32), result

  b = [2.999314767377287d0, 2.993844170297569d0, 2.982962913144534d0, 2.9667902132486006d0, 2.945503262094184d0, &
       2.9193352839727122d0, 2.8885729807284855d0, 2.853553390593274d0, 2.814660195524919d0, 2.7723195175075137d0, &
       2.7269952498697734d0, 2.6791839747726502d0, 2.6294095225512604d0, 2.5782172325201156d0, 2.526167978121472d0, &
       2.473832021878528d0, 2.4217827674798844d0, 2.3705904774487396d0, 2.3208160252273498d0, 2.2730047501302266d0, &
       2.2276804824924867d0, 2.185339804475081d0, 2.146446609406726d0, 2.1114270192715145d0, 2.080664716027288d0, &
       2.054496737905816d0, 2.033209786751399d0, 2.0170370868554657d0, 2.0061558297024313d0, 2.000685232622713d0]

  result = exp((b(ind(1)) - b(ind(4)))**2 + (b(ind(2)) - b(ind(5)))**2 + (b(ind(3)) - b(ind(6)))**2)
end function fun