From maintainers-request at octave dot org Wed Oct 19 10:02:35 2005 Subject: Re: Trouble with fsolve From: David Bateman To: "John W. Eaton" Cc: Joel Konkle-Parker , octave maintainers mailing list Date: Wed, 19 Oct 2005 16:57:02 +0200 This is a multi-part message in MIME format. --------------040803000007010706000005 Content-Type: text/plain; charset=ISO-8859-1; format=flowed Content-Transfer-Encoding: 7bit Taking this thread to maintainers, as its getting a bit more technical than a normal message to the help list... John W. Eaton wrote: >On 19-Oct-2005, David Bateman wrote: > >| If you want a completely matlab compatiable >| fsolve that calculates the function and the jacobian in a single >| function, without the mess of caching the results, then you have to >| rewrite the underlying fortran code to accept this... > >What is wrong with caching the results, provided that you know that >the function and jacobian subroutines will always be called in pairs, >in the same order (I'm not certain that the minpack routines do this, >but I think they might and it seems it would be worth checking)? > >jwe > > > The caching is quite invasive. I attach here the patch I made last year for it so you can see the type of changes I found necessary to get it to work. Note its against an octave version from September 2004 and still includes patches for odessa that has since been dropped from octave, so I don't imagine it will apply cleanly, but the changes shouldn't be enormous. It also doesn't include the changes for multiple argument functions, but those are trivial relative to this part of the patch. Also read the message at http://www.octave.org/mailing-lists/octave-maintainers/2004/627 for a description of some of the other issues with caching particular to DASRT and DASSL and issues of consistency of the interface for DASRT. If you can suggest a simpler manner to do the caching then why not... D. -- David Bateman David dot Bateman at motorola dot com Motorola Labs - Paris +33 1 69 35 48 04 (Ph) Parc Les Algorithmes, Commune de St Aubin +33 1 69 35 77 01 (Fax) 91193 Gif-Sur-Yvette FRANCE The information contained in this communication has been classified as: [x] General Business Information [ ] Motorola Internal Use Only [ ] Motorola Confidential Proprietary --------------040803000007010706000005 Content-Type: text/plain; name="changelog-handle2-2" Content-Transfer-Encoding: 7bit Content-Disposition: inline; filename="changelog-handle2-2" 2004-08-24 David Bateman * NLEqn-opt.in, DASPK-opts.in, DASRT-opts.in, DASSL-opts.in LSODE-opts.in, ODESSA-opts.in: Add "jacobian" option * ODESSA-opts.in: Add "inhomogeneity" option * DAEFunc.h (*DAERHSFuncCached, *DAEJacFuncCached): Function pointer for cached versions of function so that ires and cj args can always be passed. (DAEFunc::DAEFunc): Correctly initialize cached versions of the functions in the constructors (function_cached, jacobian_function_cached): Return cached version of user function and jacobian. (set_function_cached, set_jacobian_function_cached): Function to set the cached versions of the user function and jacobian. (fun_cached, jac_cached): The stored cached version of the user function and jacobian * ODESFunc.h (*ODES_fsub_cached, *ODES_bsub_cached, *ODES_jsub_cached): Function pointers for cached versions of functions so that the column for the inhomogeneity function can always be passed (ODESFunc::ODESFunc): Correctly initialize cached versions of the functions in the constructors (fsub_function_cached, bsub_function_cached, jsub_function_cached): Return cached version of user function inhomogeneity and jacobian. (set_fsub_function_cached, set_bsub_function_cached, set_jsub_function_cached): Function to set the cached versions of the user function, inhomogeneity and jacobian. (fsub_cached, bsub_cached jsub_cached): The stored cached version of the user function, inhomogeneity and jacobian * DASPK.cc (ddaspk_f): Pass cj argument, select version of user function to use (ddaspk_j): Select version of user function to use (user_fun_cached, user_jac_cached): Cached version of functions (do_integrate): Set up either the normal or cached version of functions to use * DASRT.cc: Include "lo-ieee.h" (user_fsub_cached, user_jsub_cached): New cached versions of user function and jacobian (dasrt_f, dasrt_j): Select version of function to use (DASRT::integrate (double)): Set up either the normal or cached version of the functions to use * ODESSA.cc (user_fsub_cached, user_bsub_cached, user_jsub_cached): Nes cached versions of the user function, inhomogeneity and jacobian. (odessa_f, odessa_j, odessa_b): Select version of function to use (ODESSA::integrate (double)): Set up either the normal or cached version of the functions to use * DASSL.cc: Include "lo-ieee.h" (user_fsub_cached, user_jsub_cached): New cached versions of user function and jacobian (dassl_f, dassl_j): Select version of function to use (DASSL::integrate (double)): Set up either the normal or cached version of the functions to use 2004-08-24 David Bateman * DLD-FUNCTIONS/fsolve.cc (fsolve_user_function_cached, fsolve_user_jacobian_cached): New versions of functions to allow user function to return jacobian as second return argument (Ffsolve): Update docs, allow function handles and inline inputs, destory temporary functions on exit. * DLD-FUNCTIONS/daspk.cc, DLD-FUNCTION/lsode.cc, DLD-FUNCTION/dasrt.cc, DLD-FUNCTION/dassl.cc, DLD-FUNCTION/odessa.cc: ditto. --------------040803000007010706000005 Content-Type: application/x-gzip; name="patch.handle2-2.gz" Content-Transfer-Encoding: base64 Content-Disposition: inline; filename="patch.handle2-2.gz" H4sICBK5NUECA3BhdGNoLmhhbmRsZTItMgDsPf1z2zayPyt/BezOy0k2rZDUt3Rpm8ZO63du kondvr7p3WgoibKpSqJCUol9Of3vbxcAQYAEJUpWnI8XTyYiQWCxWCx2F7tL8OjoiFSfTL2B P4ycd+6Tlxdnb+cn/iIKq9686gfedck2zfqJ2T4xm8Syu2an22hUzfiPHJvw/NHJyUk+mBQE y+padgbCkfqH98SuGXaD0FvsAG5rLQJXJ48I/J2/PL/q//7s4rcz8pQcdrvh2yAi5dOfLvpn ry/PL169rBzSepdnV/2zP16/wVrld86UfE/MqlkhPxC86ZK8lmcvT/uvXl+dv3r56JjAP3FN yMtnv9JOJ87QH3jO/JCWnr563j+/OvsVbl5MnWvy/saNbtyAwH9kGcLFeDkfRp4/J4EbLYN5 SJ/8N4dBnJA4JHSH/nxEnOB6OXPnURVgXcmVhlhxGvpk4BJ47gbuKG64cAMngtK3S4d2woqj wJtfk5kDv7cALPLJj0N/5H4Yh/70nbuqUsxxqBL2hFz972s6QG8esbGpxDZZoUxZoCUWSkQ7 SrHW6bPL1/+4P2spYHZnrU7TsDttibmgoGaZeex1YlXNPHZ6qvJTXPXr5p+REy7+elj2OXsB 46/e3Id1OAjR2q4T4Jl6rWs2C7FNvWY0OoxpCCnBH1A7jMhzf7qczX93h5EfPCa3Iz8y+POR vxxMXRIZhF8NJ5UePKOsxPEh5Xe+N6o8OsAyYCCYZlI2KwYB9mAXgRu6IKGiYOlWyAeyygCA ize/XNLrcQrQeGdABhYD59C7SR7YSRGwjEz89jFxUsCcKvzEAJ0q/CRAnSq90MF9THzKtH4A LKTpg3ZBPlA+gJlrtuPVfa+ZQ9aN7hbuyB0rzUn5KKHec2d4444qMVapbgwAgn2UEL98RKRK CTawph4TD2iiInasIvYrXbAMJT6FxVBiNNuEE6tVBKmt+BxK+kOKpnig3O9tHey1o+Lr5N7d 7rKO5E5pBX4Hj3hroiBC2yV19rcMdQK1ZRotycRrdYx2M16lcR+C0LH6Y3zEmfQD14b4uMdw O1baMb4XrcU484DwCgjrODVOIEA/wUJltniUnKJ0Qp6ScU8/7JbRrifDbrcNy6wntgf+cYSO ohsv7PFSeXRZjMTIsmOv0FaI3nGCXlyfYnksusXpfkqQB+VSGZljjkxCHsH93M7ob5grqKfO lSKjsmDWT1rCr3mTpsFLXbPp2cMV9JRM9LPXbhodiWkts2lY9s7Tlz/YLF2y8yitXIrvzvO4 CPwIJL076uqXXk8z173cpRbPR+78igq8sydPwFD1QrChA89BBQLXV2/AXgQreE4QX+IPJoAg PqDzHywRW4M4aO+CSRuDccZg1kLpHWolN5iDUT5yIofcgEE7vHHm1+6oCji5gfcOUBlOnTB0 Q7Bw79DI1u9aqsPhPTcsAEFpXO82OkX3Ks2GkJGUWGEEBvowZqJuV5on3Cb0+WRpqsVTR6vR +RPVgFZkPu8pPUAZ3Iyoqd8fx6KdKfjHJPJmrkHkMnKEDTOFIyS1AbqmRJQnjw2pyjRydGbE ES0kR7hC2er86ezn85f985dXZ2/e/Pb6qv8/51e/wDbi+RndNlzSEeAUINVakirZJ9WO10Dj LB7DlBZCPmi5jbo2PuHsgAn30SZIK1VrRlvsrdhfNFv0KdJVd+rOSNmrgESjBeRP71+KWD3g 1SmiUCumP1hTMRCDV8DxcvrgUJL9mAdUxBK6pVe3D+2G0UlJ+C2RQ+AppoitNx3eQglsQt9g pjZCcqcgwPJB7oEUunmzai3DqtfuNXN8u4KtFqMYZWrFrkNZ3kuPcQOGnDiB5mYPfv6OK4ZM jo9j9EUVj1XxeBUPq1BFXrcNq9G41zQrI+mpMy/Z1XyaPggLHPEaLWezu54oylCjEE9QIBJT IG4a5vhEZNYyEJq/HVnDsT/Q5VxWxsTFP8HNTxNJmlh1lZ5akdlxoqLGDkzGRlKrFAdUYssf VLNhWx1lZ5KPIgKJe4IOcLJLH5iETS9vzSDELNOxlEhm+teNR9+YkctMwRNFK4E5Z5J8bM1c nMxsh3nzU9pmcih+62dI6+AFo9y2LDl6AFPYEA7ekmIQw887Z9qLy7HD0kpYzUgVvfw+xu6P sRFleBSYlBK0iGtIuiClYsVnxFqUj9KQSfnWYB4X6l5hGhcBxa6eEqMEPKhO3fl1dAPEIgdP yW1yW2HV2G6Bj7d8NFwGgTuP+sJA7btB4Ad9sIxHUzegjUrlQ2pWdGFYaBN4YQRNSOj9G9DF 9c3kH1regIA3WoKN/Y4OKDyUEORyzb0OqEeZdSTtSfRzcBzPwTGbeYkz6WzTOSowDi1XWA3D ti3hwdUxAOtTEsxjv1yvyNLyB2KRLs6nxJaCi+ON1xj2JK6wRu1G3bCbMvdt0bHQHuQ//0mv vMq22KQ3OK9Ozy51jnV0jZvdRh12ORu3OAmMlGfdanTr9WKe9YbRrCc2RMaPiV30J+FykOvA 5OuN+29JviMzunEjJxH7MfIZ9yR0xtxyA3E14VeyFyxpT3GkzcZpMOPdwBi0uE9bDvJhDnaG yW9o40mRDiY5HTBCx/eqa5I2dKr4KyA61QG7nfBbNrca2Dr3otKP4uZvGG07MeAeiI3WhAQE 4RPJvqHfDaEB2q1SC1XPkNap9Eg6EKBDZ7B3dEgePrlxCTERWyKyYT7UamlMdl3uBpFmMKkk FySudPXBniTFx8RgByHzYOjsIJ/2jtt+RJsUbVHwY42TeMtAfSiv1Ry8WQdiHd1HhOo0cwe2 92ZLcneDDWU17C3c3WkZyAhQMCCTNOPO/WNlYDQAogWX7VUTBZEep8IgOY5yMR55ngbFxzMo MJ5B/niUXrPjGSjjGWwzHpmlJsXHMykwnkn+eJReNdENZTyTHQMZknCB/3pJGZMxahkTNbQs n3fl6IaeGzI1MvRNnLwrbot+585H3lgXgnhzlZ84Ve/W7G69UBwiAZOGUAMghex0y6yDOGjK 0S8osNr5iVNr0qbkpKn/FylTQfTAGXdvrrSxq07XanZr7WI8o8Su2sRqd2utbr1TiF1qTaMu OWHwVmiO77z5cLocueRw6p/QzXv15jBdHgI5XWeWfvJ26UVYdqzW9lzXZVUfJTYnWoHlI0r9 /ng47y+ioJIKl6hBkaP0rbwXEFUwFmIot0cVfby40RF7WxqZMo3mVpEpLqA2hqa40CoYm1Kl VMHglCq4kkZvrpK+rp7T8GyC2JAPIBvP0uZGtEVuhOyDpy5p5n2nlzrnu7LZwUZjlC6Jf43t +3J83boQSAxD8rE9lQIi3OF3YvdYkMjSBol2xjwVPpDtRG4JPnlC/vjjD/Li/A+QjniFkms5 h26mdyiExn4QBSC/Rr4bzv8WkWvvHQpGKgWhg6VL/HECajg5OCCXbkS8CKUWtPOmU/canXtY 1yAgBKMbJ6LNhYPbm9N7ilgC6z00RYnpXc99kJhS5CE7J8kWcO3UGIQ7+s7nY33E6+GmW6sa LRM0YXOvgTA+lMnaoXycYJgNo6nZew+GcT853qWiY1n21pFhI7MUiIx+QtJqvdJ2w7CB1NsE wNBKLRYBS9XcIgSGlJFiYHVLdmHfNwamWPZFg2D6RkWjYIwSZhpi0TjYepKXtqJ3HiHM/OGa 62JhyVRpGaxeM+y6/LZDC+ayY2uDYQWjYPJq3ToMpomErVMLxWJhKcm932BYED1kMOwjR8Ea tmE32pkomIhC5Qa/kLd3iX7hyzWmnYl+reswUQpJ1EtminuGvWBnc3mxdlNd75qtIhukBEwW QsMstqmuWaBra9KmGgvqjW+b6iKb6nD6kJvqi8tXp2cbGKe22RujgElDaHRNqxjj2G1gHFNi HCyo23mMY1GAWeahw/7aWQVgjB70jUf0+10+yzhgWGjcbHbtZqHIOgWRam13urVisqXZMloi OY9vv2P3bLebuEZ1roZUvYGoN1hbbyLqMZfEAdEkrsJchKGDiavoOSFz962Ryj9NXr3imavM e9LCbGjO359qPPn9xvZa1sdysAYLpdWgSKtJulXKNbMPcmuzdNvo4xBbXO1ep79wAr7dkTIZ ncCZsQ0jXPLd4lqvTWJ1K7szIwG31e6d20eTJVSL/L47H079ECRFH8gQubcRtcuRwzqw8TVr iQj9aGN8aP/OBJDawsPjzW/8mX/tzl0vuovlqQAmu302e3syW528KTXAalnr3PkUXKFXvbYB 2na/S0FymMjDnTz8ImAWaH2/qyA1vN5GB9DXxf+TvfH/5DPg/0bLsJq1B1QFg3WDNujsfqS1 0KoZVvsTaISB3iWaJspmnuK02cBVn4TAWgdFs2PYbUv1pCdeL40XRWz0eV2dszSxoJSEkazL dJCuP1hff5KuP8mtLw+CeKNxL+0zFQhQdx7UAMBWL+MAYg9Mxp92u2PUrKbqqt+JXNt5bXMI WsR3m0MrfdNBTtPB5qY6n+9gL37g9dxU2oaVStvwUQG/8WBbl/LBFnwpvHGKcCrEq/qDeOzk IJ6Ua1vm3Jkb3fgjLmEOw8gbjw+5x/sgdgAn8ZzYhc1a9cfom3hKbIv7tmPMNFV4KLdmt4xa rZYJdtwLLa0j816oalMurKZRs+VoEsfDC/1FJPrjCBKio11J7cqyekmUIMYoXSemXMMyag05 v2Gr7nNptCNKWgq1YHLbNYXnEAnYaztzsPX60OnwL9av+paeLmIhbUhumcoUieFpSbw5SCGi bXXLNuo1Ndq2I4697EJm8SKYgU1jONgc/NFAEbJYAiYZthrJWiCSwfmlfMhcGV0y95mPMVwu FlMPunvzyyXcDAJ/GXlz9+CwwoPMpXWRlwM5HHEgB1u2nje9/2+tr7jRtVsFnYAZZ7Fdw1ws s9Wt1Qo6izsGbOtkZzEUNKQ93nqXMPn6XMKMlT6mTzjPj65sLu9DOd0ulY06uvGCEQBDMrqU jpSKAgpQYehMh8spEkwBJUESpOZ0Rni5lEYgtFcS+O9h6w31VNo/JOW1Eb9sNmQbtnddy+7W 28WCfdlsyLpZNHn2s8qGDKefKBuyaTQlIjRamBz5pZ/ToT2qI+/kp9pne8qFaXTsj3XKxZeY q3iPUzk25So+1OkcYEdiev5XcDqHzZIkHi4hMXNcx24nchTIRfxMDuJoAYVb7c/5IA6r3TRs s7G/gzjWnYvx4Md1bDyIY4+HauxyuscuB3FYbdhkdKTXH2wLpjA5yPsTHcTxhZ3FEU6/urM4 rE4blnJ941kcBU/i2CIbEROZ662NZ3Hc5ySObZISw2D45PTi9OTFby+f0yPSnrDD2Xc4dXAN KAVKu1trFz1+sN4xGtLyxVtBOxC7P0/9AbDcwqenLVKO9BfIYlBId62w6fDm0qGeYgNaTYxm zjVCRh0xtOMzCoH9sj2h1ZfaF1ODT94cT3LcComVn9czbo34HoKN8xeoRt7DPyfA0TgDfwmC ZuZce3MnuGNWZ8g7x9dEA3+mYvdDMtyB7085INpRAkbiDk1dZDZNXb0uN0GXd7R+dTAsJEH1 PTHJ48dobZTNStUL+/GECb+W8GQf6JFWWg/92WLq3vZxvxmDiM9W5UIGgaDFXz5kxO4yyxuL EqALJ4g05DwR3r4MPx3Gii3uQ0NdJu4k7SIEAHNh8nEwb0WfbTvUMLCtjQncn6Tq3H7GJE0z 4a4kzT8KTk6E7lhGrS599SMjqldCPlyCznOuXSqCxC5vFK9Mf0wLxUJHhRmPNxEHikHAACAP SW6BX2O/GH2muhl0rW+5DpafwS0XM7GtoVqD2pNQ0Bpgb8anrJZYVdNnyFNxzxiTuo1fqhfz oY6Ka/P0i/fUi0EUE0Z0IYX2KW1B7QO2IFXD5P18vMMTw8FKIWVLbYsozmGNWLLZxskKW7i5 QViDUs4Whu9guG2Cmz4Ds94BWdz/8ZYygniocEjKs/gZIlemLfBBjxk3so0p4RXbjgy8uQb8 aAN4ZfyqlpEJkaUsLCCU9UqLk+9Hfn+2nEZe35uPQDb4C1K2Ddp3jEdimaalk6UzR3lFjcRC K5RWvI8aYM15j3uVWrKwUmHnKoC4Gp8XvEmWBdRST03io2JWtZBgwrjmVGEmL7PGwTpM4Emk /zv118S9XwfegouAeDB9FxdpTI94RKvsLFnFZ0mrWayHniWtbslMl1655E8X2+3yAWVUTIH5 AQD7mZ9kihLxS27VrZVOACcbrET2bNP5Kn2aS2YLx5/8Gke/ZNUj3AMFVQ8X07sqHXFm/Del 803pfFM635TON6XzsEqHC+AHUjqEfMeYhby4fHXx+1n/2U+v3lzBkP/5iAWP+QVMPVzlH9gr AqanZy9+e9k/vTiNhRkTPgaBmbv2l/TjXYQcnhydkMi9Rb8dget/zrGbHwEVZLnxnHy48J0R /crIC84XK/Lhzx9Bcn64XRmEXmDj+HoWXq/+RZ6uuDwkZVoKqyyucGuuKtjLAfnZe+fOifQc N51zZ+aydAix/Yy3o34wi7/IyMHeriorBot+2QSDll7kwRwAFwURLgDqARMdG+oXHQn9Ye4v 9Df6YwYsztMISbgc3jBXGW8o+qXZ26sqoxf7/3ycDAY/veKQ6L1/gsE26oVmyR1OEDh3bKhj LwDtET/HgYcMkLL5/hFWzc2H8YqM3HAYeAPmUKPTCWNmqFNfXQYQy5voGLXEEfPVMIUhMwjz nRM/YCwwhVXEgBVkofXcw0AVZiEd9zAQ63jogLD/layjhJtEahEXHpm8L0wjIe9vPOgQ6wI6 DF6qGrKl4uxlKE0kYlTJOf+I0NAJ3QQSr8kbrrjnGgcsE6TPisMVmS2BuQeUKp/1KtGf+oCH MFiZrFdcJUBnmDYb9TRfNGDz1jUfNRP2IIujpcpF8FTyjipG67gP0KEKt1aVLPPwvRfBTJdp nWrgvw8VdymdOKvL3JEKHu5tFDjD5DtoIlOUQjJIrLrgqt+XLNr+oSGcm4eC6GgAqNVI+bZC i+X6PeLORxk36SBwnb/irFOKsc0xFvnOjBn678SLP2hm0yE702mfPQ3jaDhvwywryZaoJJ5f Cri0kSAYp2Mm5UaClLajR0lLDoZ8Sd12sETuGOcM4+TgTeeSWVt63NF+k3BX55OBVuvSEbAH qWaaaRXP2SwkcCpSy+zCwCerZJpWbC5X2eBi1ko8ICmSydZTOlFC2TXAbpEuLUu3cWBp5Q2j 3mlk0sqLCIF4KxlGo26XyzRkCJQ6BhmzH6Qvv8Kf3tai4zhT/35ihXIgXV+w9RCedq5dwRQF eovHTM3KqRJyyojEyzQVO51ZosWawRb9JtuEA2F7Z5puaJbKEkniSTkS9CD1xoUiSVO+Ae3A SvF6jWcbQ/Bz7+3S7Yd3s4E/ZYVlVZqwVclQxra8YbIy/7yj7PIvLFRqVZ3FApYhIMX7q+gf H9J1LINJKqbnIVe2ZLVEhqe1wr60SggYT+NuhFpLpLt7UucuhywfmyRcEyY6LNaGKcW1rTos 6XWhqlR2pP5G+m+cgSJzUHQWtAq70EQoGktWJQloRqxSLK03UkpZzKVJhkxcmcrPORHKcR8V zcOYQrx5XGWzVVDSmgQp5ZPHpiX2v16fl7KqKiHbKlexSzx/IHgeNpfOchpJolbW5BYFXD5k mwHUYOGNv5yOcIMh3lmBRWFrNxJSaL6gUSGxAtV6qjd4j4aG9tUHE99Rbql5yLHlgabG96RG Ucg6/mDv9zf8jsdssYShLGBPRKAB3JA7NzpUtf3LC/oKwHzqvp3TgZZ1gfaK+kpcigVE42rO p4N1ARSVVi8vzt7OGRY0e1CAlF4LpGW0C769VJQwt9caZs1o2Jaa9PzxySaSoXOkCLclVErp CaTQvKQluKDGDkRflzW7DXpxEl9P17QgOgqM1d75Qf+yLywrW0rPa9jAMOkseRrMmo/9nn4x Qy8o6TgCzJmF1ZOPGwk24OqHe6+58344dZ1Apk2iFo/VTOwCjSV1cazmWGW1Pk+Akjcks/Ba DIVVn+GLZtf5OUe1mmXUatIXDRp102jU21J+KOtd0GdJGcb/C7cPymbJFvtyll6bLMb/CkE9 AW7VIdo0TEySVIgOc9jdECbkLpW6FPlkEbjvUANEN+4M9AOONHCB+UaitRQAUANJFfktCLFn wL8FUCwCFqbUkRz8eWmp9JuAe8lKjSGlvyvQ7NZbn29SKvvU8ifISeXfeP7SU1IbJv28bpLc 16hR4fWFJ/etez7yI6kKj+vzh/lPhhPpEf2IGStn7ynp0gkZh2yRTWhov6FG33Jgn1BL3u6k bz1kPvw9nOyYkIgiUyIPLcPfpJiWRdmvAeWmMEqH86coFQcpczJOnBAoFoFKrg6dhTPE957Z UR6IkVxWNBFFyUUp1zAfYzhRM1TKNpZGqUKL5pJwdpHKWY5JUqhoVCEXiuRyKEJku1QORe9m kzjWRI1pp4eVQm+YJO/vRGAhsAi7mgkUZ5bg488xoST+aO4Xlk9SKuVnp6XSDD5pJsg9yPtp EkG01GPMa6sUk8QWh2gDhcA6U8HliTldR6n9P5vgzElt2LbY+t2Ya6IIdSHL1IdRLPrU4uFE CEr1AR+uGOf2eSuyBFopGqFYtqSkWYslS+6gWbVKtWCqZVF1KhUgtZ/yvrbP0vymM7/pzG86 85tivJ9izE9N3LfO/Lo04wMovyRr8/TZ5et/7JC0adYN2xInEPKzoHCDzh3QPC8LM6XYZd9c iZwmFg6BRiwjStcSCCg1hru+krXIM87isArPJZPS7zzpUDDqn6GJU6oLIfq/9s61t23kCsOf 1V9BB4uihsZZiRQp0cluu+h2gS4WWKDJt7owfJE3chxLkJQmaZD/3rnPOXMhh9TNSlhgG9sa jmbOkDOHM895z1xUdXP18MClsuRGOpejEfMDbZoKuV/pjQkbmJMk3D/XAh57w/ZsFLYHEbK/ TT9esYde/sYaSyv/wrZUNMhGQP/lz2vBGHIscpyRbKLPCo7S7HQ1tShHQWoQBEKC4kLELDBC oqq4Ydr1CPnp5oKkqZZA0LSjtzI2UOZLr5aqkyvwkRifK9HCFTXM1RIN0JY4xZNtcYonkHi1 hvDd/HZ2x5YoJfsn7b6efhQ/ffed+PdXOgKfLy7YUsdA17vk4oJ++TIxf/r4RT8jk5xkRuNq LxZfh/HXEPlKP1gwuYhHTuzaj8SuDeVPIZdy0Qyz5T4ckFGqLakeETbK8DFh/9936NtE9O5m quYZ2X9lD2UGxdsyQdc/preisg+z9Rt+UR29K7x3Be8+F1fLBm3nQehv60Hoh0eyL2372lwu amM6i8pSK5QRa/H2S7JaTG9md5+ksCNnwEU9rIYZBMi5I7RSTZwtRZZaOi3+d7oyLb+bv19S u6svlDc8/5IZuxosinLRYKJV9IO1fMT4V8nDM2ZyaiPqut2S4B3HHlWV28pFFX9MMnX4Rn97 mRQurazf3jBxqI9x9sMqw1YEUWX5JSSRLhvHW83bpx9VFo4vKge1ik5lgSeMLdfZBlHLVZbp tTdMLb8Mdgxs0gzeSlH0stOHKHYZFMViVDf3rTlmXSfEwpwnZiOKGRsOvExUokWiAvtVvAZr LgYkz4JYc8VcsRnUHDfD9O3ShyWasQCNXiT9QLPd5nieGdqmAc58shWcuaJXDSBd1PldcrrW I42BUXsIwtyuu4rshWWGs/HuTRQ0z85N87SZZu8o7AZpDo5F5Gj4VvXmZDNYXlqAzXD53RHX 7F+oVfla36HncxyiEGfIN1vrfM9ZtOLo5ljG2FROf0cz1+lOvQAv2ZeXJCvGjgaeOnz733Q5 p19BX1cuxbuKvdJK6Ub+gqe6VgUSY1uzQiF81HOqeGqlsqZ2ogak3ZNywEoKmLWa8CadQo9j itFRswAqkrgYkdzsDm5mDF+3LTZY997baZsK9tlPxoE1tmKVT+G0Sp/lhgkq1cZwS3yV1EPB Wx1if56bEcmUBiS/ByYTkpcTDAerIDL51v4D/yIlFUvbpu8I2b5LUfB44WDWv0g0uMzJaADS axfDESlSTVf3eorzlSdwAcwXSRnvgffVtqlBfvXxUAXxu1x7id/x+SjzZSmvIH55TSiP9eA8 HZ1nkzjid0zyFGw/jonJXFPBzS7po/6itsh9fZGbuxfVhC77pntp+4PwuX5BzYIMS5jzh04B 2WjYMbd+5tZftcCAlmu5PGzG/ySR/M9OodpmaqDHhPzYsA9/LCNZH162FvPZEuOzXDdjfI5C Uo33qyNgdwj6tDXv0yFgt4XfqAfoSeM3+imP1ulEyyJeee73s/J05Gm3DHXLUAeVHt1a0wgq 7ZahWAb0X6/lhm08ATqakDQv1M6rIKXYxtEXnWyPU0BcyBITIuv5YsFuN4OpfLhasdTXRCAn 7DXz3dWn5HHOxA2Wy+lqwZAYngPwkwZr5MupPIKRpFaHiGblhIwGhXr7P85xmc54WtxalNQB DDdBSbWO6lNFSYsBSYtJh5I2ICT9nKSGG4eT/WGih5ZINTeYpZVqxiiCtlyu66VSbXP70T+6 fMA0QKOMzlsmDVANbPrk2M6GdyYkPrVshT2iYsK7soQmxKK9UjdCIt4Il1dsPlclV/SGma+m yXI+X6/4/SpRzOT2vZxQpzBX4HPVGmpL/W3++c2f9Log2RhJ1DBYa4BPoZxtaJ1Nqt+KzpKv UYbLkhv99h9u7gyjpV69EJr1/ffJ399Mb94mfGmgfdXI7fMk+WnJ+y8o2+Xsjzd01XqkPu3D 6q+hBNDZmGX9xklWo9CvIHOqsAcMQ2nkTJrBhxhgKob6iApMZNcA7NRLVqoyFcipj02wpNJk hCZofWpazw9tfTSNbLfD09Di7L8QePoZsGkVZvEgKmHbJGEoN8JCIVxTHhpah+OsQh9DK2/s KoJ26O2IRZ96yVNVrJ489QNUgHCStZ0acgs+o5BCEeSJYVAsAAvdNUBlTx+sq1cWrrG3icQe PEfV6aGtV0fdszv9PIImvDDXsuzzrzFRog+VTkHaUNrQRb//An+pYl/om7xsJHz8E/TA81mN FqQzyuwWnJ/LOsCtwF66T3Vq1DxLSZGOESQSNTltzqTCWQze8F6GE8zUFWgpJkDM3VZzyRc7 TffnOu7/xMf913QBMayRfJ+cUCC0Jq+uF9cFBQOwn6+EB2Nz9XbRUhOciLSqrLr14gBMQKFB uril6erM9mmb9kJ0ZMx6vAUb6YW1l6BojqQ2miO0pnqWohot+hbj0avEVCsR1Wg8tcb2vhU/ 3vwwWkQs0ycWwRuFpOo1GSC1NVSqLhIAU+3PI9hU7C/HBrZwi0VK8GqToYAU6Bw4TvyJy6Zq R6Hn8w16WH7X9gp4+zcT3/0S6x6ou0JCseIwqJGz8ELXDm4wxEciLBK6FdZqVolA2oxLE6/E WvbFqg8+t9Z9sQ7i9RV4L+pCvEz/5dSRpj3RV/FbBw4zHd/pxwWd5Ph7o7IPf5yY+KJ6lU6f 4VpVR8GftNHM+7Tf4Dd3Hj9iL/6bH5obkNEYSC4WRUmKceEESumHLwq61efRp+a01jnorK7H d9Pdw5vuT+bJoG9M9DGWYaXhvk5SMpqACOVikpGiBFK4gIOshV0N3hpEW/kbazXdGiRbfde6 cKt3rcX3yyZ4az4oSD4EcOB4mJJxqndm0PzGG90eWY3AVd8/fpg93l4ulvM1/ez8fPme3i5L sTj+os+w5J1h7ReFANaH1fx26gVYJwxgbSBZq2pClZRM93aYRQGs+ZDkUGWY/lpuRbKWfUzX rU+JaGIFycoLHELAVnzx0QvYZgUZjiBMm5UkHY+/EpjWQ8p6tWbFYDbSmlXQ0eZ06yYYa0z+ YIH2RHI96LaOgXvQBXshfPg3foWED+/XQQifnpT5OybYdN8m3hYA5No6GjjdAewDnqVY2Afz PM2wHfjoNs/yDmbpWN1Sd5ZuToJWz8/VfGc3P3fz85OYn7tJeNNJuJ7E/CbnZ4NV/vbq959b 5EMflgMyLDU9IAymdLw+rdb01VtLO5JkNRdvb0jfEUs7CrYDKkIeiJLk8ATbOUxWs9tpLSyJ gKnGsJ2QsvTQdpCVTAuSZSjG9MCmPjD4GD1COx4cv+TskKSDUgPHj3O2ZCuZOovf4/VrNO8j YvYcZG9H7KNB9TYkzMA4S1YPj6ICyzioJ7/huWDLfnV2ktBYiR5XDRh/TLIxyQyat2PLfwWc pPAxg5zkCejpxgPkZ/NLksLkWdmEjmBZeODK2R1LA7dltlK1vL/lO78fdec3UNTECppcH9Wz CEj8XXyBWAtUo6XsptjjU1exp4Mdw7JuUEvOrtgrl7ja0lMOi2LmYVHMFAnd5cFk3ur8QL/x YAU7vVO8G43MM8PKnDE/C7bCOW4+S3BWYuF3cYrNvK8xiu3MUciUsTaoID/8Xp+qz+B1zlH1 GThZPkt0s1PZ7M/q0t2oY1ZaxUYWglZx2L44o8QoY+oKXL0yeANFAQRO+6OUMUFR1YWWepi6 JgggOE9HSA8TIwgufADc/a3AB1FiXHh4QBMilLfqBDfTMckHk7DgZuw8tJn+Ztzs1bdLH1Z/ 8wTpb2p3wI8I2m2O19+EtmmgvymVkzaX4KzoWANtSdT/an1JB3lsJzO59iGOPXckwlqT9kK1 FxFOOOnvUITTkXyMWKS2YhInufqTEt/0Wn9X4pse2c2oUfC5C81lN8Gq0kJ2E67tO5PdXDti m7XuSM/ni0TCjD2ot2n5ET1nkYrT2/SozWKPgv3P9ip6G/KMTYQ+Ta/o7/igZbvuhl/Zk76y Fnl7ZU/aNhCH4SEaLGVPPKpVeJzn3A33nhsmoYW07KOt9yhDMVgW8nF7vU7YxReJpzf1gp0h u/R8RqnjBv2GaSTW6W1RtEqntwlhec6IcfLH9w1JNoagV56RvMhqRDcZuVXNIR6bzCbrUaTM JtMkhOG++Tgn+cQjsylPw56IzGaUxqY+JarQ2Fw9uIjihCGK2eR8kDfR2OQ1oUrK88HoPB/H aWyWBG4Msl/HW0EUn1fLa64eKqFE9kXbZxLF9x49k5gPyDAHhPeQTjfZsOwEPm94WlQrUXK1 9Ker+UnvkEagY+PUvzy18f71P83GSxCh5FO+N9u0ySL51cm0mRkhUqvNXLAvwbbVw7eUG3j1 8JVCnk9E4q2tfePgojaGbiMpulGC4NA0t0vROPUQP2nROD3TNKdPwdIZS582Xjo77dJuUewW xeMja7uVb4dYbXBRjFoMu8WunULqq9+ao7zpYEzSod7DPsp07ccogcr2Zyc6b/VRmn27hK+o LG6YwiN0cvCE5zuTVh2WJE33Jq16ti0By7NtCVieNRWwPFNjMbvTap8+mdVsUpCslMmI/Kwq nSNR2qejyDpP/d8jyTq/gTJpZC761cPXkIu+HFRjt1Z+aYy75V6ZF/GS5CSOFmche0tNr1sR TE0PJMmoj6TEI+UrXm1ielXu+BLTV1rGVgsLWqbX3jC9KnAVvZp7E9PrGylWvwv3ITYxvSq6 zcT0ok5LCQw/L5smpgeGA/77wRPTByeOjfPUR0w3fbv0E8tTL5bVYJ561OZGeeq1bQ6Rpz7U q2YJ2E3nD5mnHg1BZDJ2Pu3sK0+9npwPkqe+bknZkmm2n6fe0nHePFO9Ow6HyFRfPx6+Zb5V pnq13rTLVK/X40Nlqq92Jno+T6JNpnq48PecZWvrmepl5UKUE5xy7NQt8OdIYJm60y1mqrfB FTdTPbB1jdyjfZLnZKp/9VtMGvM9pKEHT1llGnqvdXo+00SoYXrM0zoFfYA3apaCPnDy6k9B Hzl2fm6TvhVPcphfviR5WXb55SPB19FgyNNxgfzyeTC/PDuvOibwVZ+phMDX+S0zj5d8Lc+H +fkoiyZfdVVWevncK/Hp10IlkKJkvxZ16eXFt1bnl5dl7iPKXFclmFffdf/iKaWXT7MBSTNg tzQbwqiFI6RPrzdTzKy4knZnfeXHU2948coE9PIGiMtAb0AZcDLtK8vbRERT+AfbUuiEf2Df wf/Ifti2dKeZ5dnONjs2xjgH+HBxtUxkK3AZtFIIOyQ/yrNtg9eIcaGmX61fCgP/qEq7tfAv U/gIgHHUPWDWRlOcmshbnrtA8BJcRt7hsLshyTzcRNw+w//YbbMaJgqGGzVwGoW+Vs9ktTAR oHGsQWF1sfdVEWYiqvOiRhlAjSwcIOr6baNK4tu+Qr0/0bFjTLm8dzHW4b7t+wTUWC0TpvEm vHYNmO7bgLPHN/N38z+mj1PGaFZa8bqxDa+9t2sar025FUwLTkztxBMTa703a51Faomp/IfE OD0tiC00kUbzycgntJyq+304VW0FdRu7U22Vdjt3qnOnOneqc6eemHbyN+UzNaO8O3eqc6cO 605V7FRdH/FOFfvA3g5X1v4BNMLnfl1bidHkdS/FUtNgMA7rwHXeWuetdd5a56113lq3w7U3 lyywmna+2ga+mn34zE1O//ru6q1ccqSD9hOLPnkpFrAf/5y8u+ThKAx4+OxnQIYlydKxIpd4 lMyU9lYHm8hzXP4dN/PlcrpasKAPFtw2Z7FquJwI5ZBEogyW+sYTksgfFiaEsciGpMg0ZvME bN5lJsGj5H1SspRkIx1DZAKpTBYAUYO4ic2dSCdvasn1dLlyRsVEk5ngtatPIhUGE/58FNFk rBaJ88JQTdjDk+TfF8/uLp69SC6e3Yt/ri+e/UcVxAF2qgIdClmR6qEqTk3eDiJYLSLNQyjD gw7GVV9hQsVEaFpkhSpnCVwLxZoP4zPZEzgakSK3whA3HMs2aU3Yc8NALhZJyDoqDepLeiJD IKNHSVTlhBQeyyj5eUX6BOIM4SUZDwaYV0xuF5dpdWBp+0dPRYLCR68ffvT6ThN8EZh4Upb5 qHQ0ZHX6kudwynuth8Kqc+HUKKdW6372Qo7phIwmAzXved/7QJpxvrshX0tNtJLC6Jy/XKP4 pKjYJL73sdSxDjJgR9LocndA7UcYVF/uu5iW+Uh4HSwpfTQeT6euAeGSnmQcsJQQsSbJwkrI 4c3GoeB3uGnBm58yl1j8mLk9aRL8gTBrHwnL4eQa2+BgihoDBdK4RBmpKs4Q78LosM/PJngG 3GneWAcT6ODpAgqaDLcf5Cvh7ceRkuG2W7tI9zxa0AywiYUCD0c41Ei8hPk7co1jP01Xbu9u F7Ar16AjJLk5VQXQ1eH4GHFDqbpMeIz9vIv4JmwFfd29uM6eNXA0iJ5g2MQ/Hg3JZGRpMofn pH7DwEsiVqdL9bOOxWw+sfXBVfct5rqNYzDxrAHuPm/QImxqRTAlDowAVgleY89z7nzWIO69 shfulFAf0QbmJRin5bcZ9lzkV5owxlBeE5LcfrxdbJzdRD2nVo6TpC66cmfpVvgM6Mm48sWZ k+/Ck7Kz7kYE6SXBzRewr9xo7C/R3nYwnhGX8JkrUAJaTIUz9m0j+WzUb2sjvEElA0v5zzx8 RsaVqh+zc/GZtp/PyZBTlONmyEvUvx5Xo2/CV1mHWg1Nr2pUelUD0osai17dMHidoPiR4HVg H6YvA3DVJmx9RCnyVsAeb0VUKSriCSz1fY6NZIJK++4iHzAV1qjQ5ooILVXX6x9svwmtE9KF 6oPgXC46apZx2kalQ1ptVHGJNmtv5tiUTePGqKKAsqj5wlPPx84krqpShav8vj4I1MVun+Os EN1Xr2VlwyzXTV6BHZo+CNnt48GA7psphj2kvozy7Seglp47H+mDX5N2jxf4/ed/vHr1k8qQ 099a3j3tGXpDyUYZyXO90cfOQZYzGeN4yd6b1amJHafMWguieFHEj/sehpxfExi6en8NXDgE uLryHMjx1nVch+u4xrG/tVG7zNkuBwUpU50rsKU9Qm2/80f7hrtwB4N9kRFxrG+1Ie0art0a aswYFyMM69AsjS8WzA0Rrr303rq03lGt7KL3O67Dwce1N4/34SoKko+HZg+vzEpS5sOjTa9j LjbTcPXl1hoRF6bMNr3eUkOw445Ls2uX/DXJkvMk9Yd35uWIx3MbY49HpDTimZbeQn+joOOq gGPZy/ePH2aPt5eL5XxNPzs/X75/vLxbigX4F3MsKu8x633+//97X4XrcAEA --------------040803000007010706000005--