# 
# L-DFP example with M=1
# 
# Written by 
# Tamara Gibson
# Applied Math Program, University of Maryland, College Park, MD 20742
# gibson@math.umd.edu
# 
# Companion to 
# Tamara Gibson, Dianne P. O'Leary and Larry Nazareth
# "BFGS with Update Skipping and Varying Memory"
# Technical Report CS-TR-3663 UMIACS-TR-96-49
# University of Maryland College Park, 1996
# 
# 1996
# 
# 
--------------------------------------------------------------------------------
> \
# Initialize Problem\
\
b := [1,1,1];\
A := array([[1,0,0],[0,2,0],[0,0,4]]);\
x0 := [0,0,0];\
H0 := array([[1,0,0],[0,1,0],[0,0,1]]);\


                                 b := [1, 1, 1]

                                     [ 1  0  0 ]
                                     [         ]
                                A := [ 0  2  0 ]
                                     [         ]
                                     [ 0  0  4 ]

                                 x0 := [0, 0, 0]

                                      [ 1  0  0 ]
                                      [         ]
                                H0 := [ 0  1  0 ]
                                      [         ]
                                      [ 0  0  1 ]
--------------------------------------------------------------------------------
> \
# Initialize Algorithm\
\
x := evalm(x0):\
H := evalm(H0):\
g := evalm(evalm(A &* x) - b) ;\


                               g := [ -1, -1, -1 ]
--------------------------------------------------------------------------------
> \
# L-DFP Iterations\
\
for k from 0 to 20 do\
  print(k);\
  # Compute Search Direction d_k\
  d := evalm(-H &* g):\
  # Line Search\
  alpha := evalm(-evalm(transpose(d) * g) / evalm(transpose(d) &* A &* d)):\
  # Compute s_k\
  s := evalm(alpha * d):\
  # Compute x_{k+1}\
  x := evalm (x + s):\
  # Compute y_k\
  y := evalm(A &* s):\
  # Compute g_{k+1}\
  g := evalm(g + y):\
  print(g);\
  # Compute H_{k+1}\
  beta := evalm(transpose(y) * H0 * y):\
  delta := evalm(transpose(y) * s):\
  H := evalm(H0 - (evalm(H0 * y * transpose(y) * H0) / beta)  + (evalm(s &* transpose(s)) / delta)):\
 od:\
\


                                        0

                               [ -4/7, -1/7, 5/7 ]

                                        1

                             [ -6/35, 9/35, -3/35 ]

                                        2

                             [ -1/35, -1/35, -3/35 ]

                                        3

                             [ -1/53, -1/53, 3/265 ]

                                        4

                                 26    27
                            [ - ----, ----, 5/4187 ]
                                4187  4187

                                        5

                              172     301        6235
                        [ - ------, -------, - ------- ]
                            154919  1394271    1394271

                                        6

                       3148675       45710075       1454905
                 [ - ----------, - -----------, - ----------- ]
                     7294051277    65646461493    65646461493

                                        7

                6440067768109      1183329214537      5383936272427
          [ - -----------------, -----------------, ----------------- ]
              36336169840374909  12112056613458303  36336169840374909

                                        8

       6891595293454210976      5656785076851084449        11559749863667868317
 [ - -----------------------, -----------------------, - -----------------------
     90598268253064400648121  90598268253064400648121    90598268253064400648121

      ]

                                        9

         4188139660870241412771899698         8665169201809534940321550287
  [ - ---------------------------------, - ---------------------------------,
      339885529515672102993133356793899    339885529515672102993133356793899

           3929999440649309123279650675
      - --------------------------------- ]
        339885529515672102993133356793899

                                       10

                    7707870486748920375138583983995670631721
             [ - ----------------------------------------------,
                 1197639617845737114271360507675902023544716363

                      2477415947626861308778694803040280644481
                 - ----------------------------------------------,
                   1197639617845737114271360507675902023544716363

                    13183355783614542642394490241246934607501
                 ---------------------------------------------- ]
                 1197639617845737114271360507675902023544716363

                                       11

             98977278561538269599561215724426982977434519509581423
      [ - -----------------------------------------------------------,
          28692253043387646854748368129098843244862795627086203894299

             2612892938614410230835505723488352847649735698579330525
          -------------------------------------------------------------,
          1061613362605342933625689620776657200059923438202189544089063

               534496969974920609901709467493613943211751497533230185
          - ------------------------------------------------------------ ]
            353871120868447644541896540258885733353307812734063181363021

                                       12

[

-

    1287541571523769318559179385116197576303030351880291329296612267514900050
 -------------------------------------------------------------------------------
 3513107514603949381339599048986337850906945724653459881092256650921890575291513

 ,

-

    2743845963976615156830459667617471850428422617267409817473411979952979815
 -------------------------------------------------------------------------------
 3513107514603949381339599048986337850906945724653459881092256650921890575291513

 ,

-

    1071454659778868808534119914559542015478070872437400771968376119242483835
 -------------------------------------------------------------------------------
 1171035838201316460446533016328779283635648574884486627030752216973963525097171

 ]

                                       13

[

-

3238007988628481633778144917870096623208939904880731188035145225263206447553866\
04468668678200
--------------------------------------------------------------------------------
1617405149750603811069568107344133138795374605781698431041738791246061846563545\
478791684956638030197

,

-

6272963331594672658252049458658989118084240012762838235256221891574620812190811\
4692042103865
--------------------------------------------------------------------------------
1797116833056226456743964563715703487550416228646331590046376434717846496181717\
19865742772959781133

,

5580555394459319705416944750334178548584056357490242477230825490708940506032063\
68220460271455
--------------------------------------------------------------------------------
1617405149750603811069568107344133138795374605781698431041738791246061846563545\
478791684956638030197

 ]

                                       14

[

-

1111799725064259006850783131718232828071488639681322517009831401220442081080294\
308357685586592129923724221760196719075
--------------------------------------------------------------------------------
1231732848502732314063045426561109812811412723674610617870801330952135648765604\
3557956152242667769059557967706891857600795789

,

1986150419169629274089222557933107566728126795846121722858580852092874989416764\
420958567652917255595151273476198459375
--------------------------------------------------------------------------------
3695198545508196942189136279683329438434238171023831853612403992856406946296813\
0673868456728003307178673903120675572802387367

    ,

3990182372629968456136871707175105916299438457791592495270236901325077834449397\
67793164816934379540702641435875445125
--------------------------------------------------------------------------------
3695198545508196942189136279683329438434238171023831853612403992856406946296813\
0673868456728003307178673903120675572802387367

 ]

                                       15

[

-

1822387079818297321532741780452067095296527327802584709542319038317210750878525\
24385326210541372238831622993464024935881871766051185193211578528825
--------------------------------------------------------------------------------
1063627635631028970865208723812422782928333993541684902004984233200453053304603\
0772444493782517265242797307771019777016748271031873328910996069686169368377

,

-

3997815433528933147241827999335153193351340948459297090340455384189221603092516\
92527483987023596973182924508150297737870084988905619111957984725275
--------------------------------------------------------------------------------
3190882906893086912595626171437268348785001980625054706014952699601359159913809\
2317333481347551795728391923313059331050244813095619986732988209058508105131

,

-

1682901667036811893793020484580132008300425628969866613880532553213642244366546\
269123382376988759476191301457899831483481989823650744069941870348475
--------------------------------------------------------------------------------
3190882906893086912595626171437268348785001980625054706014952699601359159913809\
2317333481347551795728391923313059331050244813095619986732988209058508105131

 ]

                                       16

[

-

1189116796461224751902916306659559686031997990768121469588523049652089744240831\
3917308787991098880763172748395471810962469245499829524047450292235305207322556\
00662908451226658018900
--------------------------------------------------------------------------------
2227082350836525163740440462839451996571478156738030817983449907477331570244196\
1235196700797782518501107774064838208569757893769711442046835967166936751502638\
9950445912602925818118968884877

,

-

1013624502146193200432551027631446910346941718107311292087663431088179527503306\
3684320242041088452173289946779831780092366457191583956232966746962076505383313\
47073079797684697494425
--------------------------------------------------------------------------------
7423607836121750545801468209464839988571593855793436059944833024924438567480653\
7450655669325941728337025913549460695232526312565704806822786557223122505008796\
650148637534308606039656294959

,

8695542167019856952510197930786165956874422773332557421642143170144482010087284\
8910467949865269778650881123556116510665830882162086465240174766780300008242530\
6384610996844276034375
--------------------------------------------------------------------------------
2227082350836525163740440462839451996571478156738030817983449907477331570244196\
1235196700797782518501107774064838208569757893769711442046835967166936751502638\
9950445912602925818118968884877

 ]

                                       17

[

-

3434694432944097487698254644481213720488717931744348970853492862593005893692788\
0230070231811898899935702063391999512537288030807241530024051308558493013303134\
51231755135562310206943216224633520929641554906763249642930
--------------------------------------------------------------------------------
1620099665553712043070379773216723639709471298453237751874349681158197018141882\
7661706397780270804103795535665334615667295965368864834769688442059708012073092\
02834469788705475144087289134081481729839243160098947933504287864021

,

1712003298342035683532772452792160353615635827400249207054942481576995985248889\
3834161073407614267061661336988353523731748425648813492892968594809138614270170\
23812260121060689203659821455976390183225093820165572752255
--------------------------------------------------------------------------------
1620099665553712043070379773216723639709471298453237751874349681158197018141882\
7661706397780270804103795535665334615667295965368864834769688442059708012073092\
02834469788705475144087289134081481729839243160098947933504287864021

    ,

2928615639651019487656240365898278540400199544162899387384480234770628243576033\
5772907147596026447281371637040824197457186633516837859502097118957998883200625\
12274622980971551797973169415301638784321064313976922753435
--------------------------------------------------------------------------------
1620099665553712043070379773216723639709471298453237751874349681158197018141882\
7661706397780270804103795535665334615667295965368864834769688442059708012073092\
02834469788705475144087289134081481729839243160098947933504287864021

 ]

                                       18

[

-

1023374383731228412208743455950276428087438148823820916424265585731645012196238\
1384250609154083338328918964988536288336551866317477246666382221985323008409064\
3393814149435477168235293837014102753464926093121471968516861894556438558034809\
5405368522588944729835
--------------------------------------------------------------------------------
8966864030065413269325213250191741729218593212605305334255604825993738282486731\
6294251343422611619224110982572051219947910859099872436112631446258178227546122\
2280126552639282770529131776079048693472146876086870472777636754672233446352361\
186124985465429006045637681057

,

2778286213726805151589003463479293687515935284566604080313219405017936422515278\
0117678589759698993265731186904177443768678848434805074848020267923123892909226\
8994281791967465983446343821233293034263590762372748206295344702732435603961348\
12267390131618534129515
--------------------------------------------------------------------------------
7083822583751676482766918467651475966082688637958191214061927812535053243164517\
9872458561303863179187047676231920463758849578688899224528978842543960799761436\
5601299976585033388718014103102448467842996032108627673494333036191064422618365\
33703873851768891477605376803503

    ,

-

2601673815925312767276052626297834477923700345364713142762335977039511544820074\
6991613617960911693995936258064773786036130086340083202906795411461043189567999\
6943692300008979630889904829482676552230551889583333407555559834038492405027835\
9867655640558997855495
--------------------------------------------------------------------------------
1727761605793091825065102065280847796605533814136144198551689710374403230040126\
3383526478366795897362694555178517186282646238704609566958287522571697756039374\
7707634140630495948467808317829865479961706349294787237437642203949040103077650\
0822045703701680279941594556183

 ]

                                       19

[

-

4051027107487347136376549598503681734899211268195717725306098372300702169452826\
2792309013324464175029114978249891758242396747065090117782798801255708216626765\
4714871526226501887099226420814693187248894530857305789545130427547174667948281\
507437847447036012281844697794288901437185844264868750691676592670625
--------------------------------------------------------------------------------
2520902386579953514949096934304352609841934727711335388334653706214939988210681\
6178108580879581626316571359514561373384771818883863185055110190518305902341679\
9125432021493184209748495803616201934408033098840556308521105711204915716359721\
1602084604762601205525305423102135682206972183819497141113608174804027805866997

,

-

3473487469681021817982298827203731445199501587393259764438041259937367196283855\
3720832114194617659280103710618103031744824583534733709215380847252034633430137\
6999874827335071110486256334743084765431794696045559642860788889963710884555253\
892062541178505066086604537348806371731440526834508252057265852357275
--------------------------------------------------------------------------------
8403007955266511716496989781014508699473115759037784627782179020716466627368938\
7260361936265272087721904531715204577949239396279543950183700635061019674472266\
3751440071643947365828319345387339781360110329468521028403685704016385721199070\
534028201587533735175101807700711894068990727939832380371202724934675935288999

,

-

1028522166378367263049337573124253374998245493633398971372667391079640172896208\
4453817032506615939198749930982451554466474910863502894408770112377113793481625\
7587364256306621411369239253222831556910511577862851903066574750992733214055389\
8015285121265246442306839580839583314369007348243123904331254382795
--------------------------------------------------------------------------------
4272715909457548330422198193736190864138872419849720997177379163076169471543528\
1657811154033189197146731117821290463364020032006547771279847780539501529392677\
8178698341513871541946603056976613448149208642102637811052721544415111383660544\
33933637368857647551276363103426028511982579386771137984976409742441149251983

 ]

                                       20

[

-

8860817059540554932815768531637134267675095422837930811231744769796003752104343\
9682988458340845354071936793578464565958604224130945233798190950179084656612286\
6771060984695292190215313227897303901039580932889731654482630081172531466700438\
6989961817205976264665159187659478566941757466223343664710815250503001373595033\
645053013251387805479638018330130536039550
--------------------------------------------------------------------------------
1608152843713164805181650286930841891703119654494303214794725043346921135681097\
5552931712512270161853007283607271597520540929307432949158090231493657207552028\
3808631690105475272746293217600710127678216747603036797019544839690918667054228\
4291877810696902946683308457529992596087478098620145889120685978925104513554311\
63707241356050945554767175781917187087349404709414909

,

-

2084371507060195857761255623002687770685227793439772493260895964873876860290028\
5773532119356283159372409336435270763998748057982123440160824127429411830134839\
5484133898571140352172480303462610617528229687182833063786647936307023838443319\
8371318157129868091742288423765859471602733370403442315604300147422163084710998\
47028482826895374474876706270272874449475
--------------------------------------------------------------------------------
1608152843713164805181650286930841891703119654494303214794725043346921135681097\
5552931712512270161853007283607271597520540929307432949158090231493657207552028\
3808631690105475272746293217600710127678216747603036797019544839690918667054228\
4291877810696902946683308457529992596087478098620145889120685978925104513554311\
63707241356050945554767175781917187087349404709414909

,

8603380847363693009709084604412200126100238130696578777049014008835661191284219\
0767033666633199739199638238200467887056321640143585696527491173087222414071857\
5943964551527637454820988722954431237521058115577907386527466596294287457275721\
3631625479725267840882915452929388760864986000432527930319973549422070630421844\
32352944215764017211212369106847410575625
--------------------------------------------------------------------------------
5956121643382091871043149210854969969270813535164085980721203864247856058078139\
0936784120415815414270397346693598509335336775212714626511445301828360027970475\
4846784037427686195356641546669296769178580546677914063035351258114513581682327\
5155103002581122024752994287148120726249918883778318107854392514537424124275228\
285453383557442427954339843774710632864792767015367

 ]
