Sherman-Morrison/tests/qmckl_test_fortran/qmckl_test_f.f90

program QMCkl_SMW_test
    use qmckl
    implicit none

    integer :: i, j, col
    integer (qmckl_context) :: context
    integer (qmckl_exit_code) :: rc
    integer(kind=8) :: dim, n_updates
    integer(kind=8), dimension(:),   allocatable :: Updates_index
    real(c_double) :: breakdown
    real(c_double), dimension(:,:), allocatable :: S, S_inv, Updates, U

    context = qmckl_context_create()
    rc = QMCKL_FAILURE
    breakdown = 1e-3
    dim = 21
    n_updates = 3
    allocate( &
              Updates_index(n_updates), &
              S(dim,dim), &
              S_inv(dim,dim), &
              Updates(dim,n_updates), &
              U(dim,n_updates) &
              )

    Updates_index = (/ 12, 13, 21 /)
    Updates = reshape((/  0.102125488221645E+00, 0.962643101811409E-01,-0.180085301399231E-01,-0.116052091121674E+00, 0.100100971758366E-01, 0.168905649334192E-01,-0.168554261326790E-01, 0.334977582097054E-01, 0.775983557105064E-01,-0.509594380855560E-01,-0.110948169603944E-01,-0.301160980015993E-01,-0.210931494832039E+00,-0.409047538414598E-02, 0.239971160888672E+00,-0.458180941641331E-01, 0.231322161853313E-01,-0.548814050853252E-01, 0.185034181922674E-01,-0.714888703078032E-02,-0.105082295835018E+00,-0.727039668709040E-02,-0.329451821744442E-01, 0.230633586645126E+00, 0.322856456041336E-01, 0.188164815306664E+00, 0.976034477353096E-01, 0.124376058578491E+00,-0.241902604699135E+00,-0.214029267430305E+00,-0.164923388510942E-01,-0.796175077557564E-01, 0.552074126899242E-01, 0.794003605842590E-01,-0.999749153852463E-01, 0.135955372825265E-01,-0.876737236976624E-01, 0.210409909486771E+00, 0.177621953189373E-01,-0.151124328374863E+00, 0.247122272849083E+00,-0.161797225475311E+00,-0.817385390400887E-01, 0.108162150718272E-02, 0.693925749510527E-02, 0.493063451722264E-02,-0.334013439714909E-01, 0.323923602700233E-01, 0.120032556355000E+00,-0.119067849591374E-01,-0.334841385483742E-01, 0.125436363741755E-01,-0.145120620727539E+00,-0.101416520774364E-01,-0.723919421434402E-01, 0.134314149618149E+00,-0.125767393037677E-01, 0.361425074515864E-03,-0.351854860782623E-01,-0.506495498120785E-01,-0.291761625558138E-01,-0.932136957999319E-03,-0.500183776021004E-01 /), shape(Updates))
    S = reshape((/ -0.123240642249584E+00, 0.150952935218811E+00, 0.871514901518822E-01,-0.150895461440086E+00,-0.871202647686005E-01,-0.123183816671371E+00, 0.100359566509724E+00,-0.135056734085083E+00,-0.744080543518066E-01,-0.120210982859135E+00,-0.757721215486526E-01, 0.525853671133518E-01, 0.102125488221645E+00,-0.307190138846636E-01, 0.202103536576033E-01,-0.211704343557358E+00,-0.152789223939180E-01,-0.396802611649036E-01, 0.156803593039513E+00, 0.273103892803192E+00,-0.666790679097176E-01,-0.383377403020859E+00, 0.469519078731537E+00,-0.271040737628937E+00, 0.469401627779007E+00,-0.271042704582214E+00, 0.383267551660538E+00, 0.116775326430798E+00,-0.936452969908714E-01, 0.369134694337845E-01,-0.204805508255959E-01,-0.173394829034805E-01,-0.197960838675499E+00, 0.962643101811409E-01, 0.274471193552017E+00, 0.145023554563522E+00,-0.179225616157055E-01, 0.136946782469749E+00,-0.174913182854652E+00,-0.128193218261003E-01, 0.225567314773798E-01,-0.171542761381716E-02,-0.109934650361538E+00, 0.393143796827644E-03, 0.155300989747047E+00,-0.378105294657871E-03, 0.154906913638115E+00, 0.109414055943489E+00, 0.162110984325409E+00,-0.104241691529751E+00,-0.199010908603668E+00,-0.178463548421860E+00, 0.115249522030354E+00,-0.751317143440247E-01,-0.180085301399231E-01, 0.776397660374641E-01,-0.162087708711624E+00, 0.115154005587101E+00,-0.111836232244968E+00, 0.215133726596832E+00,-0.165389522910118E+00,-0.160914268344641E-01,-0.184740740805864E-01,-0.505007863044739E+00,-0.618503987789154E+00, 0.357068866491318E+00, 0.618444144725800E+00,-0.357080936431885E+00,-0.504969716072083E+00, 0.778413712978363E-01, 0.487910546362400E-01,-0.115782171487808E-01,-0.530001707375050E-01, 0.214109313674271E-02,-0.194653034210205E+00,-0.116052091121674E+00,-0.260605067014694E+00,-0.137633457779884E+00, 0.422291792929173E-01, 0.133351176977158E+00,-0.170262515544891E+00,-0.304439514875412E-01, 0.533868409693241E-01, 0.596583448350430E-02,-0.947399914264679E+00, 0.440122239524499E-04,-0.133992493152618E+01, 0.392828042095061E-04, 0.134024262428284E+01,-0.947782099246979E+00,-0.155572161078453E+00,-0.790303274989128E-01,-0.178044617176056E+00, 0.141150355339050E+00,-0.857104435563087E-01, 0.108973577618599E+00, 0.100100971758366E-01, 0.567631311714649E-01,-0.203381881117821E+00, 0.391849316656590E-01, 0.140034869313240E+00,-0.172185719013214E+00, 0.583750754594803E-01,-0.916111283004284E-03,-0.292119309306145E-01,-0.351941259577870E-02,-0.429157400503755E-02, 0.200041988864541E-02, 0.343890837393701E-02,-0.242624944075942E-02,-0.279381335712969E-02, 0.133628234267235E+00, 0.164535552263260E+00,-0.303669814020395E-01, 0.335082374513149E-01,-0.115301951766014E+00,-0.643811970949173E-01, 0.168905649334192E-01,-0.131173312664032E+00,-0.228150542825460E-01,-0.140278100967407E+00, 0.121666260063648E+00,-0.169033091515303E-01, 0.515771955251694E-01,-0.180635631084442E+00, 0.195857465267181E+00,-0.520476046949625E-02, 0.580822164192796E-02,-0.399843230843544E-02, 0.576473958790302E-02,-0.271833199076355E-02, 0.426918268203735E-02, 0.147665038704872E+00,-0.117447517812252E+00, 0.147420719265938E+00, 0.135893806815147E+00, 0.422651544213295E-01,-0.695310905575752E-01,-0.168554261326790E-01, 0.921243280172348E-01, 0.117380328476429E+00, 0.138013571500778E+00,-0.623580329120159E-01,-0.136343047022820E+00, 0.131682857871056E+00,-0.132132634520531E+00, 0.119739077985287E+00,-0.241280291229486E-01, 0.295925457030535E-01, 0.161907169967890E-01,-0.280807800590992E-01,-0.170978084206581E-01,-0.229465700685978E-01, 0.170763313770294E+00,-0.246933296322823E+00,-0.118787530809641E-01,-0.233876928687096E-01,-0.235260799527168E+00,-0.413223467767239E-01, 0.334977582097054E-01, 0.146476894617081E+00, 0.599658722057939E-02,-0.962516665458679E-01, 0.239439934492111E+00,-0.289920177310705E-01, 0.406377874314785E-01, 0.131046742200851E+00, 0.136250451207161E+00,-0.887189060449600E-02,-0.451377010904253E-03,-0.123341614380479E-01,-0.489434518385679E-03, 0.119013069197536E-01,-0.825028307735920E-02, 0.120959408581257E+00, 0.899301841855049E-01, 0.183146566152573E+00,-0.16
    S_inv = reshape((/  0.311919238733319E+02, 0.391211616261786E+02, 0.230468185418855E+02, 0.385544395586345E+02, 0.233591113858616E+02,-0.312005871898339E+02,-0.751758654676229E+00,-0.440356424219148E+01,-0.157135293238612E+01, 0.311608268850865E+01, 0.269147143852639E+01,-0.295675853033276E+01, 0.695722472661512E+01,-0.656548036745546E+01, 0.357345834272637E+01,-0.127226560851769E+00, 0.211510373691347E+01, 0.225402683839516E+01, 0.138095041843308E+00, 0.144755926374198E+00,-0.129650521353452E+01, 0.748881862107261E+02, 0.913799202443389E+02, 0.511079244308798E+02, 0.908925159786654E+02, 0.524296869618370E+02,-0.743781724291279E+02,-0.308858192791382E+01,-0.883914805747968E+01,-0.251520306199258E+01, 0.120501905090525E+02, 0.968749784123027E+01,-0.128360153025071E+02, 0.124602848636645E+02,-0.108315097729706E+02, 0.538637730088244E+01,-0.155954216668878E+01, 0.733071759611877E+01, 0.106785970043980E+02,-0.526404570077434E+00, 0.126274276019119E+01,-0.176385334707397E+01, 0.981074335673308E+01, 0.142283742204400E+02, 0.106061203769129E+02, 0.142205894966935E+02, 0.106658385090178E+02,-0.999705856360834E+01,-0.235530921594611E+01,-0.397822961500586E+01,-0.615755307417915E+00, 0.302954990045915E+01, 0.260774312202084E+01,-0.324474908269315E+01, 0.506349143477028E+01,-0.403331084954255E+01, 0.216159755322281E+01,-0.356236831434911E+00, 0.229287231445069E+01, 0.345423911992016E+01,-0.522840315767084E-01,-0.798660202084652E+00,-0.669484792656382E-01,-0.509168796576331E+00,-0.626061168716399E+00, 0.106282405724139E+00, 0.187208027678628E+00,-0.361241794759462E+00,-0.152171018401898E+00, 0.786005474778943E-02, 0.335536466772946E-01, 0.204150929795949E-01,-0.431873457736125E-01,-0.266642681570363E-01, 0.354230568719380E-01,-0.662672792520103E-01, 0.541727654659900E-01,-0.215053141567270E-01,-0.313441408257335E-02,-0.360061392608368E-01,-0.364432600501005E-01,-0.281329679904227E-02, 0.147042058872862E-01, 0.983501095240267E-02,-0.515300150381945E+00,-0.414148002449120E+00,-0.483927352922599E+00,-0.414236968594294E+00, 0.892894843202583E-02, 0.163447351211741E+00,-0.263796606276840E-01,-0.422167712445766E-01,-0.244456092073667E-01, 0.343720228179847E-01, 0.188162844929908E-01,-0.227730288498809E-01, 0.389866683213947E-01,-0.313257678215315E-01, 0.166784616908835E-01,-0.748798889931418E-02, 0.141492354958890E-01, 0.166134715418269E-01,-0.569910669664292E-02,-0.139449114927599E-01,-0.277054057678700E-02, 0.320422300739978E+02, 0.396796238151428E+02, 0.225744499599300E+02, 0.385319291283625E+02, 0.228620731935291E+02,-0.315787444812651E+02,-0.201688119637744E+01,-0.691868674131645E+01,-0.337527478391271E+01, 0.737932786002596E+01, 0.532046905755991E+01,-0.704646020919806E+01, 0.116798006516197E+02,-0.103900446585121E+02, 0.367492551215626E+01, 0.541292498306359E+00, 0.688933313780047E+01, 0.643920432296654E+01, 0.640069090766210E+00,-0.275193635264416E+01,-0.188561818645725E+01, 0.236539170987535E+02, 0.295003969224203E+02, 0.174888381981832E+02, 0.294531380129501E+02, 0.175044078873038E+02,-0.241830822724836E+02, 0.314282931518621E+01, 0.403784215759159E+01, 0.169155742879292E+01,-0.554066860884136E+01,-0.525713150347206E+01, 0.394727380633284E+01,-0.641137465954821E+01, 0.423388578775427E+01, 0.306527904979264E-01, 0.202300380704869E+01,-0.459569791294739E+01,-0.468425326686841E+01, 0.121902516832242E+01,-0.707896318054816E+00, 0.441694177169657E+00,-0.964925469634583E+02,-0.118149667679186E+03,-0.672201839913842E+02,-0.117411734713197E+03,-0.686625184338945E+02, 0.961995663618476E+02, 0.451976509920424E+01, 0.100355852033043E+02, 0.362241589185164E+01,-0.130069484167442E+02,-0.110760346601338E+02, 0.139329830045770E+02,-0.153873537547839E+02, 0.146103537661863E+02,-0.765612032914714E+01, 0.886407137879243E+00,-0.705337777438639E+01,-0.116989846863909E+02, 0.592924574459365E-01,-0.301336277255331E+00, 0.221249720574218E+01, 0.408696880275784E+02, 0.498549177991145E+02, 0.283099915956156E+02, 0.498508177346691E+02, 0.288351506046625E+02,-0.408181390328127E+02, 0.274581588439313E+01, 0.442598366847707E+01, 0.258341971052896E+01,-
    S = transpose(S)
    ! Updates_index = (/ 2, 4 /)
    ! Updates = reshape((/ 0, -1, 0, 1, -1, -1, 0, -1 /), shape(Updates))
    ! S = reshape((/ 1,0,1,-1,0,1,1,0,-1,0,-1,0,1,1,1,1 /), shape(S))
    ! S = transpose(S)
    ! S_inv = reshape((/ 1,1,-1,-1,-1,0,1,0,1,2,-2,-1,1,1,-1,0 /), shape(S_inv))

    !! Write current S and S_inv to file for check in Octave
    open(unit = 2000, file = "Slater_old.dat")
    open(unit = 3000, file = "Slater_old_inv.dat")
    do i=1,dim
        do j=1,dim
            write(2000,"(E23.15, 1X)", advance="no") S(i,j)
            write(3000,"(E23.15, 1X)", advance="no") S_inv(i,j)
        end do
        write(2000,*)
        write(3000,*)
    end do
    close(2000)
    close(3000)

    !! Transform replacement updates in 'Updates'
    !! into additive updates in 'U' to compute the inverse.dat"
    !! Apply the (replacement) updates to S
    do j=1,n_updates
        do i=1,dim
            col = Updates_index(j)
            U(i,j) = Updates(i,j) - S(i,col)
            S(i,col) = Updates(i,j)
        end do
    end do

    !! Update S_inv; S_inv needs to be transposed first
    S_inv = transpose(S_inv)
    rc = qmckl_sherman_morrison_smw32s(context, dim, n_updates, U, &
                Updates_index, breakdown, S_inv);
    !! S_inv_t needs to be transposed back before it
    !! can be multiplied with S to test unity
    S_inv = transpose(S_inv)

    !! Write new S and S_inv to file for check in Octave
    open(unit = 4000, file = "Slater.dat")
    open(unit = 5000, file = "Slater_inv.dat")
    do i=1,dim
        do j=1,dim
            write(4000,"(E23.15, 1X)", advance="no") S(i,j)
            write(5000,"(E23.15, 1X)", advance="no") S_inv(i,j)
        end do
        write(4000,*)
        write(5000,*)
    end do
    close(4000)
    close(5000)

    deallocate(S, S_inv, Updates, U, Updates_index)
end program