diff options
| -rw-r--r-- | COPYING | 674 | ||||
| -rw-r--r-- | makefile | 32 | ||||
| -rw-r--r-- | src/bifurchaser.cpp | 359 | ||||
| -rw-r--r-- | src/domain_check.cpp | 274 | ||||
| -rw-r--r-- | src/domain_eigen.cpp | 187 | ||||
| -rw-r--r-- | src/domain_eigen.h | 28 | ||||
| -rw-r--r-- | src/domain_energy.cpp | 2103 | ||||
| -rw-r--r-- | src/domain_energy.h | 41 | ||||
| -rw-r--r-- | src/domain_improve.cpp | 145 | ||||
| -rw-r--r-- | src/domain_increase.cpp | 120 | ||||
| -rw-r--r-- | src/domain_minimize.cpp | 300 | ||||
| -rw-r--r-- | src/domain_minimize.h | 24 | ||||
| -rw-r--r-- | src/domain_newton.h | 224 | ||||
| -rw-r--r-- | src/eigenvalues.cpp | 156 | ||||
| -rw-r--r-- | src/evolve.cpp | 237 | ||||
| -rw-r--r-- | src/gradget.cpp | 107 | ||||
| -rw-r--r-- | src/hessget.cpp | 107 | ||||
| -rw-r--r-- | src/initialize.cpp | 201 | ||||
| -rw-r--r-- | src/perturb.cpp | 113 |
19 files changed, 5432 insertions, 0 deletions
@@ -0,0 +1,674 @@ + GNU GENERAL PUBLIC LICENSE + Version 3, 29 June 2007 + + Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/> + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + Preamble + + The GNU General Public License is a free, copyleft license for +software and other kinds of works. + + The licenses for most software and other practical works are designed +to take away your freedom to share and change the works. By contrast, +the GNU General Public License is intended to guarantee your freedom to +share and change all versions of a program--to make sure it remains free +software for all its users. We, the Free Software Foundation, use the +GNU General Public License for most of our software; it applies also to +any other work released this way by its authors. You can apply it to +your programs, too. + + When we speak of free software, we are referring to freedom, not +price. Our General Public Licenses are designed to make sure that you +have the freedom to distribute copies of free software (and charge for +them if you wish), that you receive source code or can get it if you +want it, that you can change the software or use pieces of it in new +free programs, and that you know you can do these things. + + To protect your rights, we need to prevent others from denying you +these rights or asking you to surrender the rights. Therefore, you have +certain responsibilities if you distribute copies of the software, or if +you modify it: responsibilities to respect the freedom of others. + + For example, if you distribute copies of such a program, whether +gratis or for a fee, you must pass on to the recipients the same +freedoms that you received. You must make sure that they, too, receive +or can get the source code. And you must show them these terms so they +know their rights. + + Developers that use the GNU GPL protect your rights with two steps: +(1) assert copyright on the software, and (2) offer you this License +giving you legal permission to copy, distribute and/or modify it. + + For the developers' and authors' protection, the GPL clearly explains +that there is no warranty for this free software. For both users' and +authors' sake, the GPL requires that modified versions be marked as +changed, so that their problems will not be attributed erroneously to +authors of previous versions. + + Some devices are designed to deny users access to install or run +modified versions of the software inside them, although the manufacturer +can do so. This is fundamentally incompatible with the aim of +protecting users' freedom to change the software. The systematic +pattern of such abuse occurs in the area of products for individuals to +use, which is precisely where it is most unacceptable. Therefore, we +have designed this version of the GPL to prohibit the practice for those +products. If such problems arise substantially in other domains, we +stand ready to extend this provision to those domains in future versions +of the GPL, as needed to protect the freedom of users. + + Finally, every program is threatened constantly by software patents. +States should not allow patents to restrict development and use of +software on general-purpose computers, but in those that do, we wish to +avoid the special danger that patents applied to a free program could +make it effectively proprietary. To prevent this, the GPL assures that +patents cannot be used to render the program non-free. + + The precise terms and conditions for copying, distribution and +modification follow. + + TERMS AND CONDITIONS + + 0. Definitions. + + "This License" refers to version 3 of the GNU General Public License. + + "Copyright" also means copyright-like laws that apply to other kinds of +works, such as semiconductor masks. + + "The Program" refers to any copyrightable work licensed under this +License. Each licensee is addressed as "you". "Licensees" and +"recipients" may be individuals or organizations. + + To "modify" a work means to copy from or adapt all or part of the work +in a fashion requiring copyright permission, other than the making of an +exact copy. The resulting work is called a "modified version" of the +earlier work or a work "based on" the earlier work. + + A "covered work" means either the unmodified Program or a work based +on the Program. + + To "propagate" a work means to do anything with it that, without +permission, would make you directly or secondarily liable for +infringement under applicable copyright law, except executing it on a +computer or modifying a private copy. Propagation includes copying, +distribution (with or without modification), making available to the +public, and in some countries other activities as well. + + To "convey" a work means any kind of propagation that enables other +parties to make or receive copies. Mere interaction with a user through +a computer network, with no transfer of a copy, is not conveying. + + An interactive user interface displays "Appropriate Legal Notices" +to the extent that it includes a convenient and prominently visible +feature that (1) displays an appropriate copyright notice, and (2) +tells the user that there is no warranty for the work (except to the +extent that warranties are provided), that licensees may convey the +work under this License, and how to view a copy of this License. If +the interface presents a list of user commands or options, such as a +menu, a prominent item in the list meets this criterion. + + 1. Source Code. + + The "source code" for a work means the preferred form of the work +for making modifications to it. "Object code" means any non-source +form of a work. + + A "Standard Interface" means an interface that either is an official +standard defined by a recognized standards body, or, in the case of +interfaces specified for a particular programming language, one that +is widely used among developers working in that language. + + The "System Libraries" of an executable work include anything, other +than the work as a whole, that (a) is included in the normal form of +packaging a Major Component, but which is not part of that Major +Component, and (b) serves only to enable use of the work with that +Major Component, or to implement a Standard Interface for which an +implementation is available to the public in source code form. A +"Major Component", in this context, means a major essential component +(kernel, window system, and so on) of the specific operating system +(if any) on which the executable work runs, or a compiler used to +produce the work, or an object code interpreter used to run it. + + The "Corresponding Source" for a work in object code form means all +the source code needed to generate, install, and (for an executable +work) run the object code and to modify the work, including scripts to +control those activities. However, it does not include the work's +System Libraries, or general-purpose tools or generally available free +programs which are used unmodified in performing those activities but +which are not part of the work. For example, Corresponding Source +includes interface definition files associated with source files for +the work, and the source code for shared libraries and dynamically +linked subprograms that the work is specifically designed to require, +such as by intimate data communication or control flow between those +subprograms and other parts of the work. + + The Corresponding Source need not include anything that users +can regenerate automatically from other parts of the Corresponding +Source. + + The Corresponding Source for a work in source code form is that +same work. + + 2. Basic Permissions. + + All rights granted under this License are granted for the term of +copyright on the Program, and are irrevocable provided the stated +conditions are met. This License explicitly affirms your unlimited +permission to run the unmodified Program. The output from running a +covered work is covered by this License only if the output, given its +content, constitutes a covered work. This License acknowledges your +rights of fair use or other equivalent, as provided by copyright law. + + You may make, run and propagate covered works that you do not +convey, without conditions so long as your license otherwise remains +in force. You may convey covered works to others for the sole purpose +of having them make modifications exclusively for you, or provide you +with facilities for running those works, provided that you comply with +the terms of this License in conveying all material for which you do +not control copyright. Those thus making or running the covered works +for you must do so exclusively on your behalf, under your direction +and control, on terms that prohibit them from making any copies of +your copyrighted material outside their relationship with you. + + Conveying under any other circumstances is permitted solely under +the conditions stated below. Sublicensing is not allowed; section 10 +makes it unnecessary. + + 3. Protecting Users' Legal Rights From Anti-Circumvention Law. + + No covered work shall be deemed part of an effective technological +measure under any applicable law fulfilling obligations under article +11 of the WIPO copyright treaty adopted on 20 December 1996, or +similar laws prohibiting or restricting circumvention of such +measures. + + When you convey a covered work, you waive any legal power to forbid +circumvention of technological measures to the extent such circumvention +is effected by exercising rights under this License with respect to +the covered work, and you disclaim any intention to limit operation or +modification of the work as a means of enforcing, against the work's +users, your or third parties' legal rights to forbid circumvention of +technological measures. + + 4. Conveying Verbatim Copies. + + You may convey verbatim copies of the Program's source code as you +receive it, in any medium, provided that you conspicuously and +appropriately publish on each copy an appropriate copyright notice; +keep intact all notices stating that this License and any +non-permissive terms added in accord with section 7 apply to the code; +keep intact all notices of the absence of any warranty; and give all +recipients a copy of this License along with the Program. + + You may charge any price or no price for each copy that you convey, +and you may offer support or warranty protection for a fee. + + 5. Conveying Modified Source Versions. + + You may convey a work based on the Program, or the modifications to +produce it from the Program, in the form of source code under the +terms of section 4, provided that you also meet all of these conditions: + + a) The work must carry prominent notices stating that you modified + it, and giving a relevant date. + + b) The work must carry prominent notices stating that it is + released under this License and any conditions added under section + 7. This requirement modifies the requirement in section 4 to + "keep intact all notices". + + c) You must license the entire work, as a whole, under this + License to anyone who comes into possession of a copy. This + License will therefore apply, along with any applicable section 7 + additional terms, to the whole of the work, and all its parts, + regardless of how they are packaged. This License gives no + permission to license the work in any other way, but it does not + invalidate such permission if you have separately received it. + + d) If the work has interactive user interfaces, each must display + Appropriate Legal Notices; however, if the Program has interactive + interfaces that do not display Appropriate Legal Notices, your + work need not make them do so. + + A compilation of a covered work with other separate and independent +works, which are not by their nature extensions of the covered work, +and which are not combined with it such as to form a larger program, +in or on a volume of a storage or distribution medium, is called an +"aggregate" if the compilation and its resulting copyright are not +used to limit the access or legal rights of the compilation's users +beyond what the individual works permit. Inclusion of a covered work +in an aggregate does not cause this License to apply to the other +parts of the aggregate. + + 6. Conveying Non-Source Forms. + + You may convey a covered work in object code form under the terms +of sections 4 and 5, provided that you also convey the +machine-readable Corresponding Source under the terms of this License, +in one of these ways: + + a) Convey the object code in, or embodied in, a physical product + (including a physical distribution medium), accompanied by the + Corresponding Source fixed on a durable physical medium + customarily used for software interchange. + + b) Convey the object code in, or embodied in, a physical product + (including a physical distribution medium), accompanied by a + written offer, valid for at least three years and valid for as + long as you offer spare parts or customer support for that product + model, to give anyone who possesses the object code either (1) a + copy of the Corresponding Source for all the software in the + product that is covered by this License, on a durable physical + medium customarily used for software interchange, for a price no + more than your reasonable cost of physically performing this + conveying of source, or (2) access to copy the + Corresponding Source from a network server at no charge. + + c) Convey individual copies of the object code with a copy of the + written offer to provide the Corresponding Source. This + alternative is allowed only occasionally and noncommercially, and + only if you received the object code with such an offer, in accord + with subsection 6b. + + d) Convey the object code by offering access from a designated + place (gratis or for a charge), and offer equivalent access to the + Corresponding Source in the same way through the same place at no + further charge. You need not require recipients to copy the + Corresponding Source along with the object code. If the place to + copy the object code is a network server, the Corresponding Source + may be on a different server (operated by you or a third party) + that supports equivalent copying facilities, provided you maintain + clear directions next to the object code saying where to find the + Corresponding Source. Regardless of what server hosts the + Corresponding Source, you remain obligated to ensure that it is + available for as long as needed to satisfy these requirements. + + e) Convey the object code using peer-to-peer transmission, provided + you inform other peers where the object code and Corresponding + Source of the work are being offered to the general public at no + charge under subsection 6d. + + A separable portion of the object code, whose source code is excluded +from the Corresponding Source as a System Library, need not be +included in conveying the object code work. + + A "User Product" is either (1) a "consumer product", which means any +tangible personal property which is normally used for personal, family, +or household purposes, or (2) anything designed or sold for incorporation +into a dwelling. In determining whether a product is a consumer product, +doubtful cases shall be resolved in favor of coverage. For a particular +product received by a particular user, "normally used" refers to a +typical or common use of that class of product, regardless of the status +of the particular user or of the way in which the particular user +actually uses, or expects or is expected to use, the product. A product +is a consumer product regardless of whether the product has substantial +commercial, industrial or non-consumer uses, unless such uses represent +the only significant mode of use of the product. + + "Installation Information" for a User Product means any methods, +procedures, authorization keys, or other information required to install +and execute modified versions of a covered work in that User Product from +a modified version of its Corresponding Source. The information must +suffice to ensure that the continued functioning of the modified object +code is in no case prevented or interfered with solely because +modification has been made. + + If you convey an object code work under this section in, or with, or +specifically for use in, a User Product, and the conveying occurs as +part of a transaction in which the right of possession and use of the +User Product is transferred to the recipient in perpetuity or for a +fixed term (regardless of how the transaction is characterized), the +Corresponding Source conveyed under this section must be accompanied +by the Installation Information. But this requirement does not apply +if neither you nor any third party retains the ability to install +modified object code on the User Product (for example, the work has +been installed in ROM). + + The requirement to provide Installation Information does not include a +requirement to continue to provide support service, warranty, or updates +for a work that has been modified or installed by the recipient, or for +the User Product in which it has been modified or installed. Access to a +network may be denied when the modification itself materially and +adversely affects the operation of the network or violates the rules and +protocols for communication across the network. + + Corresponding Source conveyed, and Installation Information provided, +in accord with this section must be in a format that is publicly +documented (and with an implementation available to the public in +source code form), and must require no special password or key for +unpacking, reading or copying. + + 7. Additional Terms. + + "Additional permissions" are terms that supplement the terms of this +License by making exceptions from one or more of its conditions. +Additional permissions that are applicable to the entire Program shall +be treated as though they were included in this License, to the extent +that they are valid under applicable law. If additional permissions +apply only to part of the Program, that part may be used separately +under those permissions, but the entire Program remains governed by +this License without regard to the additional permissions. + + When you convey a copy of a covered work, you may at your option +remove any additional permissions from that copy, or from any part of +it. (Additional permissions may be written to require their own +removal in certain cases when you modify the work.) You may place +additional permissions on material, added by you to a covered work, +for which you have or can give appropriate copyright permission. + + Notwithstanding any other provision of this License, for material you +add to a covered work, you may (if authorized by the copyright holders of +that material) supplement the terms of this License with terms: + + a) Disclaiming warranty or limiting liability differently from the + terms of sections 15 and 16 of this License; or + + b) Requiring preservation of specified reasonable legal notices or + author attributions in that material or in the Appropriate Legal + Notices displayed by works containing it; or + + c) Prohibiting misrepresentation of the origin of that material, or + requiring that modified versions of such material be marked in + reasonable ways as different from the original version; or + + d) Limiting the use for publicity purposes of names of licensors or + authors of the material; or + + e) Declining to grant rights under trademark law for use of some + trade names, trademarks, or service marks; or + + f) Requiring indemnification of licensors and authors of that + material by anyone who conveys the material (or modified versions of + it) with contractual assumptions of liability to the recipient, for + any liability that these contractual assumptions directly impose on + those licensors and authors. + + All other non-permissive additional terms are considered "further +restrictions" within the meaning of section 10. If the Program as you +received it, or any part of it, contains a notice stating that it is +governed by this License along with a term that is a further +restriction, you may remove that term. If a license document contains +a further restriction but permits relicensing or conveying under this +License, you may add to a covered work material governed by the terms +of that license document, provided that the further restriction does +not survive such relicensing or conveying. + + If you add terms to a covered work in accord with this section, you +must place, in the relevant source files, a statement of the +additional terms that apply to those files, or a notice indicating +where to find the applicable terms. + + Additional terms, permissive or non-permissive, may be stated in the +form of a separately written license, or stated as exceptions; +the above requirements apply either way. + + 8. Termination. + + You may not propagate or modify a covered work except as expressly +provided under this License. Any attempt otherwise to propagate or +modify it is void, and will automatically terminate your rights under +this License (including any patent licenses granted under the third +paragraph of section 11). + + However, if you cease all violation of this License, then your +license from a particular copyright holder is reinstated (a) +provisionally, unless and until the copyright holder explicitly and +finally terminates your license, and (b) permanently, if the copyright +holder fails to notify you of the violation by some reasonable means +prior to 60 days after the cessation. + + Moreover, your license from a particular copyright holder is +reinstated permanently if the copyright holder notifies you of the +violation by some reasonable means, this is the first time you have +received notice of violation of this License (for any work) from that +copyright holder, and you cure the violation prior to 30 days after +your receipt of the notice. + + Termination of your rights under this section does not terminate the +licenses of parties who have received copies or rights from you under +this License. If your rights have been terminated and not permanently +reinstated, you do not qualify to receive new licenses for the same +material under section 10. + + 9. Acceptance Not Required for Having Copies. + + You are not required to accept this License in order to receive or +run a copy of the Program. Ancillary propagation of a covered work +occurring solely as a consequence of using peer-to-peer transmission +to receive a copy likewise does not require acceptance. However, +nothing other than this License grants you permission to propagate or +modify any covered work. These actions infringe copyright if you do +not accept this License. Therefore, by modifying or propagating a +covered work, you indicate your acceptance of this License to do so. + + 10. Automatic Licensing of Downstream Recipients. + + Each time you convey a covered work, the recipient automatically +receives a license from the original licensors, to run, modify and +propagate that work, subject to this License. You are not responsible +for enforcing compliance by third parties with this License. + + An "entity transaction" is a transaction transferring control of an +organization, or substantially all assets of one, or subdividing an +organization, or merging organizations. If propagation of a covered +work results from an entity transaction, each party to that +transaction who receives a copy of the work also receives whatever +licenses to the work the party's predecessor in interest had or could +give under the previous paragraph, plus a right to possession of the +Corresponding Source of the work from the predecessor in interest, if +the predecessor has it or can get it with reasonable efforts. + + You may not impose any further restrictions on the exercise of the +rights granted or affirmed under this License. For example, you may +not impose a license fee, royalty, or other charge for exercise of +rights granted under this License, and you may not initiate litigation +(including a cross-claim or counterclaim in a lawsuit) alleging that +any patent claim is infringed by making, using, selling, offering for +sale, or importing the Program or any portion of it. + + 11. Patents. + + A "contributor" is a copyright holder who authorizes use under this +License of the Program or a work on which the Program is based. The +work thus licensed is called the contributor's "contributor version". + + A contributor's "essential patent claims" are all patent claims +owned or controlled by the contributor, whether already acquired or +hereafter acquired, that would be infringed by some manner, permitted +by this License, of making, using, or selling its contributor version, +but do not include claims that would be infringed only as a +consequence of further modification of the contributor version. For +purposes of this definition, "control" includes the right to grant +patent sublicenses in a manner consistent with the requirements of +this License. + + Each contributor grants you a non-exclusive, worldwide, royalty-free +patent license under the contributor's essential patent claims, to +make, use, sell, offer for sale, import and otherwise run, modify and +propagate the contents of its contributor version. + + In the following three paragraphs, a "patent license" is any express +agreement or commitment, however denominated, not to enforce a patent +(such as an express permission to practice a patent or covenant not to +sue for patent infringement). To "grant" such a patent license to a +party means to make such an agreement or commitment not to enforce a +patent against the party. + + If you convey a covered work, knowingly relying on a patent license, +and the Corresponding Source of the work is not available for anyone +to copy, free of charge and under the terms of this License, through a +publicly available network server or other readily accessible means, +then you must either (1) cause the Corresponding Source to be so +available, or (2) arrange to deprive yourself of the benefit of the +patent license for this particular work, or (3) arrange, in a manner +consistent with the requirements of this License, to extend the patent +license to downstream recipients. "Knowingly relying" means you have +actual knowledge that, but for the patent license, your conveying the +covered work in a country, or your recipient's use of the covered work +in a country, would infringe one or more identifiable patents in that +country that you have reason to believe are valid. + + If, pursuant to or in connection with a single transaction or +arrangement, you convey, or propagate by procuring conveyance of, a +covered work, and grant a patent license to some of the parties +receiving the covered work authorizing them to use, propagate, modify +or convey a specific copy of the covered work, then the patent license +you grant is automatically extended to all recipients of the covered +work and works based on it. + + A patent license is "discriminatory" if it does not include within +the scope of its coverage, prohibits the exercise of, or is +conditioned on the non-exercise of one or more of the rights that are +specifically granted under this License. You may not convey a covered +work if you are a party to an arrangement with a third party that is +in the business of distributing software, under which you make payment +to the third party based on the extent of your activity of conveying +the work, and under which the third party grants, to any of the +parties who would receive the covered work from you, a discriminatory +patent license (a) in connection with copies of the covered work +conveyed by you (or copies made from those copies), or (b) primarily +for and in connection with specific products or compilations that +contain the covered work, unless you entered into that arrangement, +or that patent license was granted, prior to 28 March 2007. + + Nothing in this License shall be construed as excluding or limiting +any implied license or other defenses to infringement that may +otherwise be available to you under applicable patent law. + + 12. No Surrender of Others' Freedom. + + If conditions are imposed on you (whether by court order, agreement or +otherwise) that contradict the conditions of this License, they do not +excuse you from the conditions of this License. If you cannot convey a +covered work so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you may +not convey it at all. For example, if you agree to terms that obligate you +to collect a royalty for further conveying from those to whom you convey +the Program, the only way you could satisfy both those terms and this +License would be to refrain entirely from conveying the Program. + + 13. Use with the GNU Affero General Public License. + + Notwithstanding any other provision of this License, you have +permission to link or combine any covered work with a work licensed +under version 3 of the GNU Affero General Public License into a single +combined work, and to convey the resulting work. The terms of this +License will continue to apply to the part which is the covered work, +but the special requirements of the GNU Affero General Public License, +section 13, concerning interaction through a network will apply to the +combination as such. + + 14. Revised Versions of this License. + + The Free Software Foundation may publish revised and/or new versions of +the GNU General Public License from time to time. Such new versions will +be similar in spirit to the present version, but may differ in detail to +address new problems or concerns. + + Each version is given a distinguishing version number. If the +Program specifies that a certain numbered version of the GNU General +Public License "or any later version" applies to it, you have the +option of following the terms and conditions either of that numbered +version or of any later version published by the Free Software +Foundation. If the Program does not specify a version number of the +GNU General Public License, you may choose any version ever published +by the Free Software Foundation. + + If the Program specifies that a proxy can decide which future +versions of the GNU General Public License can be used, that proxy's +public statement of acceptance of a version permanently authorizes you +to choose that version for the Program. + + Later license versions may give you additional or different +permissions. However, no additional obligations are imposed on any +author or copyright holder as a result of your choosing to follow a +later version. + + 15. Disclaimer of Warranty. + + THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY +APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT +HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY +OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, +THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM +IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF +ALL NECESSARY SERVICING, REPAIR OR CORRECTION. + + 16. Limitation of Liability. + + IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING +WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS +THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY +GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE +USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF +DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD +PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), +EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF +SUCH DAMAGES. + + 17. Interpretation of Sections 15 and 16. + + If the disclaimer of warranty and limitation of liability provided +above cannot be given local legal effect according to their terms, +reviewing courts shall apply local law that most closely approximates +an absolute waiver of all civil liability in connection with the +Program, unless a warranty or assumption of liability accompanies a +copy of the Program in return for a fee. + + END OF TERMS AND CONDITIONS + + How to Apply These Terms to Your New Programs + + If you develop a new program, and you want it to be of the greatest +possible use to the public, the best way to achieve this is to make it +free software which everyone can redistribute and change under these terms. + + To do so, attach the following notices to the program. It is safest +to attach them to the start of each source file to most effectively +state the exclusion of warranty; and each file should have at least +the "copyright" line and a pointer to where the full notice is found. + + <one line to give the program's name and a brief idea of what it does.> + Copyright (C) <year> <name of author> + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see <http://www.gnu.org/licenses/>. + +Also add information on how to contact you by electronic and paper mail. + + If the program does terminal interaction, make it output a short +notice like this when it starts in an interactive mode: + + <program> Copyright (C) <year> <name of author> + This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'. + This is free software, and you are welcome to redistribute it + under certain conditions; type `show c' for details. + +The hypothetical commands `show w' and `show c' should show the appropriate +parts of the General Public License. Of course, your program's commands +might be different; for a GUI interface, you would use an "about box". + + You should also get your employer (if you work as a programmer) or school, +if any, to sign a "copyright disclaimer" for the program, if necessary. +For more information on this, and how to apply and follow the GNU GPL, see +<http://www.gnu.org/licenses/>. + + The GNU General Public License does not permit incorporating your program +into proprietary programs. If your program is a subroutine library, you +may consider it more useful to permit linking proprietary applications with +the library. If this is what you want to do, use the GNU Lesser General +Public License instead of this License. But first, please read +<http://www.gnu.org/philosophy/why-not-lgpl.html>. diff --git a/makefile b/makefile new file mode 100644 index 0000000..a098946 --- /dev/null +++ b/makefile @@ -0,0 +1,32 @@ +# domain_tools - modulated domain toolkit +# See COPYING file for copyright and license details. + +CC = g++ +CFLAGS = -Ofast -std=c++11 -fopenmp +LDFLAGS = -lgslcblas -lgsl -latlas + +OBJ = domain_energy domain_minimize domain_eigen +BIN = bifurchaser initialize domain_increase domain_improve hessget gradget evolve eigenvalues + +all: opt ${OBJ:%=obj/%.o} ${BIN:%=bin/%} + +opt: + @echo build options: + @echo "CC = ${CC}" + @echo "CFLAGS = ${CFLAGS}" + @echo "LDFLAGS = ${LDFLAGS}" + +obj/%.o: src/%.cpp + @echo CC -c -o $@ + @${CC} ${CFLAGS} ${LDFLAGS} -c -o $@ $< + +bin/%: src/%.cpp ${OBJ:%=obj/%.o} + @echo CC -o $@ + @${CC} ${COBJ:%=obj/%.o} ${OBJ:%=obj/%.o} ${CFLAGS} ${LDFLAGS} -o $@ $< + +clean: + @echo cleaning: + @echo rm -f ${OBJ:%=obj/%.o} ${BIN:%=bin/%} + @rm -f ${OBJ:%=obj/%.o} ${BIN:%=bin/%} + +.PHONY: all clean diff --git a/src/bifurchaser.cpp b/src/bifurchaser.cpp new file mode 100644 index 0000000..1ef595b --- /dev/null +++ b/src/bifurchaser.cpp @@ -0,0 +1,359 @@ +/* bifurcation_chaser.cpp + * + * Copyright (C) 2013 Jaron Kent-Dobias + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +/* A program which facilitates automated mapping of bifurcation points in the + * energy of a system where the Hessian is available. Currently, only a one + * dimensional parameter space is supported. + */ + +#include "domain_energy.h" +#include "domain_minimize.h" +#include "domain_eigen.h" + +#include <unistd.h> +#include <stdio.h> +#include <iostream> +#include <stdlib.h> +#include <math.h> +#include <string> + +// GSL includes. +#include <gsl/gsl_sf.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_complex.h> +#include <gsl/gsl_complex_math.h> +#include <gsl/gsl_vector.h> +#include <gsl/gsl_permutation.h> +#include <gsl/gsl_permute_vector.h> +#include <gsl/gsl_blas.h> +#include <gsl/gsl_linalg.h> +#include <gsl/gsl_sort_vector.h> + + +void geteigenvalues(gsl_vector *eigenvalues, unsigned n, const gsl_vector *z, double c) { + gsl_matrix *hess; + hess = gsl_matrix_alloc(3 * n + 3, 3 * n + 3); + + domain_energy_fixedHessian(hess, n, z, c); + + domain_eigen_values(eigenvalues, 3 * n + 3, 2 * n, hess); + gsl_matrix_free(hess); +} + + +int domain_eigen_perturb(gsl_vector *z, unsigned k, unsigned n, + unsigned eigen_num, double a0,double a_fact, double c, double eps, + double g, double N, double energy_thres) { + + printf("Beginning perturbation.\n"); + + double a, temp_eigenval, eigenval; + unsigned kk; + + gsl_vector *temp_z, *eigenvalues, *eigenvector; + gsl_permutation *eigenorder; + gsl_matrix *hess; + + eigenvalues = gsl_vector_alloc(3 * n + 3); + eigenvector = gsl_vector_alloc(3 * n + 3); + temp_z = gsl_vector_alloc(3 * n + 3); + hess = gsl_matrix_alloc(3 * n + 3, 3 * n + 3); + + eigenorder = gsl_permutation_alloc(3 * n + 3); + + domain_energy_fixedHessian(hess, n, z, c); + + domain_eigen_values(eigenvalues, 3 * n + 3, 2 * n, hess); + domain_eigen_sort(eigenorder, 3 * n + 3, eigen_num, eigenvalues); + + kk = gsl_permutation_get(eigenorder, k); + + eigenval = gsl_vector_get(eigenvalues, kk); + + printf("Getting eigenvector.\n"); + domain_energy_fixedHessian(hess, n, z, c); + domain_eigen_vector(eigenvector, 3 * n + 3, 2 * n, kk, hess); + + a = a0; + int failed = 0; + + printf("Starting loop.\n"); + while (true) { + gsl_vector_memcpy(temp_z, z); + gsl_blas_daxpy(a, eigenvector, temp_z); + failed = domain_minimize_fixed(z, n, c, eps, N, 0.9, g, 0.9); + + if (failed) { + printf("Relaxation failed, reducing perturb size.\n"); + a *= 0.1; + } else { + + domain_energy_fixedHessian(hess, n, z, c); + + domain_eigen_values(eigenvalues, 3 * n + 3, 2 * n, hess); + domain_eigen_sort(eigenorder, 3 * n + 3, eigen_num, eigenvalues); + + kk = gsl_permutation_get(eigenorder, k); + + temp_eigenval = gsl_vector_get(eigenvalues, kk); + + printf("BIFUR: Perturbing %i, %e, %e\n", k, eigenval, temp_eigenval); + + if (GSL_SIGN(temp_eigenval) != GSL_SIGN(eigenval)) { + gsl_vector_memcpy(z, temp_z); + break; + } + + a *= a_fact; + } + } + + gsl_vector_free(eigenvector); + gsl_vector_free(temp_z); + gsl_vector_free(eigenvalues); + + gsl_permutation_free(eigenorder); + + return 0; +} + + +bool bifur_consent() { + printf(" (y/n): "); + char in; + in = getchar(); + getchar(); + if (in == 'y') return true; + else return false; +} + + +// Initializes the program. +int main(int argc, char *argv[]) { + + int opt, min_fails, eigen_follow, eigen_num, examining; + unsigned n, N, j, a, last_pert, ii, old_ii; + double c, dc0, dc, g0, g, eigen_thres, approach_thres, eps, state, old_state; + char *filename, str[19], in; + bool subcrit, reset; + + // Setting default values. + eps = 0; + eigen_thres = 1e-13; + approach_thres = 1e-6; + eigen_follow = -1; + examining = -1; + eigen_num = 25; + last_pert = 0; + subcrit = false; + reset = false; + dc = 0; + + j = 0; + + gsl_vector *z, *old_z, *eigenvalues, *eigenstate, *old_eigenstate, *eigenchanges; + gsl_permutation *eigenorder, *old_eigenorder; + + while ((opt = getopt(argc, argv, "n:c:d:g:h:i:N:p:m:j:e:t:s")) != -1) { + switch (opt) { + case 'n': + n = atoi(optarg); + break; + case 'N': + N = atoi(optarg); + break; + case 'j': + j = atoi(optarg); + break; + case 'c': + c = atof(optarg); + break; + case 'd': + dc0 = atof(optarg); + break; + case 'h': + dc = atof(optarg); + break; + case 'g': + g0 = atof(optarg); + break; + case 'i': + filename = optarg; + break; + case 'm': + eigen_follow = atof(optarg); + break; + case 'e': + eps = atof(optarg); + break; + case 's': + subcrit = true; + break; + case 't': + approach_thres = atof(optarg); + break; + default: + exit(EXIT_FAILURE); + } + } + + z = gsl_vector_alloc(3 * n + 3); + old_z = gsl_vector_alloc(3 * n + 3); + eigenvalues = gsl_vector_alloc(3 * n + 3); + eigenorder = gsl_permutation_alloc(3 * n + 3); + old_eigenorder = gsl_permutation_alloc(3 * n + 3); + old_eigenstate = gsl_vector_alloc(3 * n + 3); + eigenstate = gsl_vector_alloc(3 * n + 3); + + FILE *f = fopen(filename, "r"); + gsl_vector_fscanf(f, z); + fclose(f); + + g = g0; + if (dc == 0) dc = dc0; + + min_fails = domain_minimize_fixed(z, n, c, eps, N, 0.9, 1, 0.9); + + if (min_fails) { + printf("BIFUR: Initial relaxation failed, exiting.\n"); + return 1; + } + + geteigenvalues(eigenvalues, n, z, c); + domain_eigen_state(old_eigenstate, eigenvalues, n, eigen_thres); + domain_eigen_sort(old_eigenorder, 3 * n + 3, eigen_num, eigenvalues); + + while (true) { + j += 1; + c += dc; + reset = false; + + gsl_vector_memcpy(old_z, z); + + printf("BIFUR: Step %05d, starting with c = %f.\n", j, c); + + min_fails = domain_minimize_fixed(z, n, c, eps, N, 0.9, 1, 0.9); + + if (min_fails) { + printf("BIFUR: Newton's method failed to converge, reducing step size.\n"); + c -= dc; + j -= 1; + last_pert = 0; + gsl_vector_memcpy(z, old_z); + dc *= 0.1; + reset = true; + } else { + + geteigenvalues(eigenvalues, n, z, c); + domain_eigen_sort(eigenorder, 3 * n + 3, eigen_num, eigenvalues); + domain_eigen_state(eigenstate, eigenvalues, n, eigen_thres); + + if (eigen_follow > -1) examining = eigen_follow; + + for (unsigned i = 0; i < eigen_num; i++) { + ii = gsl_permutation_get(eigenorder, i); + old_ii = gsl_permutation_get(old_eigenorder, i); + + state = gsl_vector_get(eigenstate, ii); + old_state = gsl_vector_get(old_eigenstate, old_ii); + + if (state != old_state) { + if (i == examining) { + c -= dc; + gsl_vector_memcpy(z, old_z); + gsl_vector_memcpy(eigenstate, old_eigenstate); + gsl_permutation_memcpy(eigenorder, old_eigenorder); + j -= 1; + dc *= 0.1; + reset = true; + last_pert = 0; + } if (examining == -1 && state != 0 && old_state != 0) { + printf("BIFUR: Eigenvalue %i changed sign past threshold to %e. Examine?", i, + gsl_vector_get(eigenvalues, ii)); + if (bifur_consent()) { + + examining = i; + c -= dc; + gsl_vector_memcpy(z, old_z); + gsl_vector_memcpy(eigenstate, old_eigenstate); + gsl_permutation_memcpy(eigenorder, old_eigenorder); + j -= 1; + dc *= 0.1; + reset = true; + last_pert = 0; + + break; + } + } + } + } + + if (!reset && examining > -1 && fabs(dc) < approach_thres) { + + if (!subcrit) { + c += GSL_SIGN(dc) * approach_thres; + domain_minimize_fixed(z, n, c, eps, N, 0.9, 1, 0.9); + } + + printf("BIFUR: Perturbing at c = %.8f.\n", c); + + domain_eigen_perturb(z, examining, n, eigen_num, 1, 1.1, c, eps, g, N, 0); + geteigenvalues(eigenvalues, n, z, c); + domain_eigen_sort(eigenorder, 3 * n + 3, eigen_num, eigenvalues); + domain_eigen_state(eigenstate, eigenvalues, n, eigen_thres); + + if (subcrit) dc = - GSL_SIGN(dc) * approach_thres; + else dc = GSL_SIGN(dc) * approach_thres; + + examining = -1; + last_pert = 0; + } + + if (!reset) { + + gsl_vector_memcpy(old_eigenstate, eigenstate); + gsl_permutation_memcpy(old_eigenorder, eigenorder); + + if (last_pert > 10 && fabs(dc) < fabs(dc0)) { + last_pert = 0; + dc = GSL_SIGN(dc) * fmin(fabs(dc) * 10, fabs(dc0)); + } + + last_pert += 1; + + double energy = domain_energy_fixedEnergy(n, z, c); + + sprintf(str, "output/out-%05d.dat", j); + FILE *fout = fopen(str, "w"); + fprintf(fout, "%.10e\t%.10e\n", c, energy); + for (unsigned i = 0; i < eigen_num; i++) { + ii = gsl_permutation_get(eigenorder, i); + + fprintf(fout, "%.10e\t", gsl_vector_get(eigenvalues, ii)); + } + fprintf(fout, "\n"); + gsl_vector_fprintf(fout, z, "%.10e"); + fclose(fout); + } + } + } + + gsl_vector_free(z); + +} + diff --git a/src/domain_check.cpp b/src/domain_check.cpp new file mode 100644 index 0000000..5270663 --- /dev/null +++ b/src/domain_check.cpp @@ -0,0 +1,274 @@ +/* domain_improve.cpp + * + * Copyright (C) 2013 Jaron Kent-Dobias + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +/* A program which facilitates automated mapping of bifurcation points in the + * energy of a system where the Hessian is available. Currently, only a one + * dimensional parameter space is supported. + */ + +#include "domain_energy.h" +#include "domain_minimize.h" + +#include <unistd.h> +#include <stdio.h> +#include <iostream> +#include <stdlib.h> +#include <math.h> +#include <string> + +// GSL includes. +#include <gsl/gsl_sf.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_complex.h> +#include <gsl/gsl_complex_math.h> +#include <gsl/gsl_vector.h> +#include <gsl/gsl_permutation.h> +#include <gsl/gsl_permute_vector.h> +#include <gsl/gsl_blas.h> +#include <gsl/gsl_linalg.h> +#include <gsl/gsl_sort_vector.h> + +void bifur_eigenvalues(gsl_vector *eigenvalues, unsigned n, + const gsl_vector *z, double c) { + + double eigenvalue; + + gsl_vector *beta; + gsl_vector_complex *alpha; + gsl_matrix *hess, *modI; + gsl_eigen_gen_workspace *w; + + alpha = gsl_vector_complex_alloc(3 * n + 3); + beta = gsl_vector_alloc(3 * n + 3); + hess = gsl_matrix_alloc(3 * n + 3, 3 * n + 3); + modI = gsl_matrix_alloc(3 * n + 3, 3 * n + 3); + w = gsl_eigen_gen_alloc(3 * n + 3); + + gsl_matrix_set_zero(modI); + for (unsigned i = 0; i < 2 * n; i++) gsl_matrix_set(modI, i, i, 1); + + domain_energy_hessian(hess, n, z, c); + + gsl_eigen_gen(hess, modI, alpha, beta, w); + + for (unsigned i = 0; i < 3 * n + 3; i++) { + eigenvalue = gsl_vector_complex_get(alpha, i).dat[0] / gsl_vector_get(beta, i); + gsl_vector_set(eigenvalues, i, eigenvalue); + } + + gsl_vector_free(beta); + gsl_vector_complex_free(alpha); + gsl_matrix_free(modI); + gsl_matrix_free(hess); + gsl_eigen_gen_free(w); +} + +void bifur_trueEigenvalues(gsl_vector *eigenvalues, unsigned n, + const gsl_vector *z, double c) { + + double eigenvalue; + + gsl_vector *beta; + gsl_vector_complex *alpha; + gsl_matrix *hess, *modI; + gsl_eigen_gen_workspace *w; + + alpha = gsl_vector_complex_alloc(3 * n + 4); + beta = gsl_vector_alloc(3 * n + 4); + hess = gsl_matrix_alloc(3 * n + 4, 3 * n + 4); + modI = gsl_matrix_alloc(3 * n + 4, 3 * n + 4); + w = gsl_eigen_gen_alloc(3 * n + 4); + + gsl_matrix_set_zero(modI); + for (unsigned i = 0; i < 2 * n + 1; i++) gsl_matrix_set(modI, i, i, 1); + + domain_energy_truehessian(hess, n, z, c); + + gsl_eigen_gen(hess, modI, alpha, beta, w); + + for (unsigned i = 0; i < 3 * n + 4; i++) { + eigenvalue = gsl_vector_complex_get(alpha, i).dat[0] / gsl_vector_get(beta, i); + gsl_vector_set(eigenvalues, i, eigenvalue); + } + + gsl_vector_free(beta); + gsl_vector_complex_free(alpha); + gsl_matrix_free(modI); + gsl_matrix_free(hess); + gsl_eigen_gen_free(w); +} + + +void bifur_eigensort(gsl_permutation *eigenorder, unsigned n, unsigned eigen_num, + const gsl_vector *eigenvalues) { + + unsigned ii; + + gsl_vector *abs_eigenvalues; + + abs_eigenvalues = gsl_vector_alloc(3 * n + 3); + + for (unsigned i = 0; i < 3 * n + 3; i++) { + gsl_vector_set(abs_eigenvalues, i, fabs(gsl_vector_get(eigenvalues, i))); + } + + gsl_sort_vector_index(eigenorder, abs_eigenvalues); + + gsl_vector_memcpy(abs_eigenvalues, eigenvalues); + + for (unsigned i = eigen_num; i < 3 * n + 3; i++) { + ii = gsl_permutation_get(eigenorder, i); + + gsl_vector_set(abs_eigenvalues, ii, INFINITY); + } + + gsl_sort_vector_index(eigenorder, abs_eigenvalues); + + gsl_vector_free(abs_eigenvalues); +} + +void bifur_trueEigensort(gsl_permutation *eigenorder, unsigned n, unsigned eigen_num, + const gsl_vector *eigenvalues) { + + unsigned ii; + + gsl_vector *abs_eigenvalues; + + abs_eigenvalues = gsl_vector_alloc(3 * n + 4); + + for (unsigned i = 0; i < 3 * n + 4; i++) { + gsl_vector_set(abs_eigenvalues, i, fabs(gsl_vector_get(eigenvalues, i))); + } + + gsl_sort_vector_index(eigenorder, abs_eigenvalues); + + gsl_vector_memcpy(abs_eigenvalues, eigenvalues); + + for (unsigned i = eigen_num; i < 3 * n + 4; i++) { + ii = gsl_permutation_get(eigenorder, i); + + gsl_vector_set(abs_eigenvalues, ii, INFINITY); + } + + gsl_sort_vector_index(eigenorder, abs_eigenvalues); + + gsl_vector_free(abs_eigenvalues); +} + + + +// Initializes the program. +int main(int argc, char *argv[]) { + + int opt, min_fails; + unsigned n, N, num; + double c, g0, g, eps, energy; + char *filename; + bool eigenpres = true; + + // Setting default values. + eps = 0; + num = 25; + + gsl_vector *z, *old_z, *eigenvalues, *trueEigenvalues; + gsl_permutation *eigenorder, *trueEigenorder; + + while ((opt = getopt(argc, argv, "n:c:d:g:h:i:N:p:m:j:e:t:s")) != -1) { + switch (opt) { + case 'n': + n = atoi(optarg); + break; + case 'N': + N = atoi(optarg); + break; + case 'g': + g0 = atof(optarg); + break; + case 'i': + filename = optarg; + break; + case 'e': + eps = atof(optarg); + break; + default: + exit(EXIT_FAILURE); + } + } + + z = gsl_vector_alloc(3 * n + 3); + old_z = gsl_vector_alloc(3 * n + 3); + eigenvalues = gsl_vector_alloc(3 * n + 3); + trueEigenvalues = gsl_vector_alloc(3 * n + 4); + eigenorder = gsl_permutation_alloc(3 * n + 3); + trueEigenorder = gsl_permutation_alloc(3 * n + 4); + + g = g0; + + char ch; + double throwaway; + + FILE *f = fopen(filename, "r+"); + while (ch != '\n') ch = fgetc(f); + ch = 'a'; + while (ch != '\n' && ch != '\t') ch = fgetc(f); + if (ch == '\n') eigenpres = false; + + rewind(f); + + fscanf(f, "%le\t", &c); + fscanf(f, "%le\n", &energy); + + if (eigenpres) { + ch = 'a'; + while (ch != '\n') ch = fgetc(f); + } + gsl_vector_fscanf(f, z); + fclose(f); + + min_fails = domain_minimize(z, n, c, eps, g, N, 4, 2, 0.9); + + if (min_fails) { + printf("BIFUR: Initial relaxation failed, exiting.\n"); + return 1; + } + + bifur_eigenvalues(eigenvalues, n, z, c); + bifur_eigensort(eigenorder, n, num, eigenvalues); + bifur_trueEigenvalues(trueEigenvalues, n, z, c); + bifur_trueEigensort(trueEigenorder, n, num, trueEigenvalues); + + energy = domain_energy_energy(n, z, c); + unsigned ii; + + FILE *newf = fopen(filename, "w"); + fprintf(newf, "%.12le\t%.12le\n", c, energy); + for (unsigned i = 0; i < num; i++) { + ii = gsl_permutation_get(eigenorder, i); + fprintf(newf, "%.12le\t", gsl_vector_get(eigenvalues, ii)); + } + fprintf(newf, "\n"); + for (unsigned i = 0; i < num; i++) { + ii = gsl_permutation_get(trueEigenorder, i); + fprintf(newf, "%.12le\t", gsl_vector_get(trueEigenvalues, ii)); + } + fprintf(newf, "\n"); + for (unsigned i = 0; i < 3 * n + 3; i++) { + fprintf(newf, "%.12le\t", gsl_vector_get(z, i)); + } + fclose(newf); +} diff --git a/src/domain_eigen.cpp b/src/domain_eigen.cpp new file mode 100644 index 0000000..c02240e --- /dev/null +++ b/src/domain_eigen.cpp @@ -0,0 +1,187 @@ +/* domain_eigen.cpp + * + * Copyright (C) 2013 Jaron Kent-Dobias + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +/* A set of utilities for find the generalized eigenvalues and eigenvectors of + * modulated domains. + */ + +// GSL includes. +#include <gsl/gsl_sf.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_complex.h> +#include <gsl/gsl_complex_math.h> +#include <gsl/gsl_vector.h> +#include <gsl/gsl_permutation.h> +#include <gsl/gsl_permute_vector.h> +#include <gsl/gsl_blas.h> +#include <gsl/gsl_linalg.h> +#include <gsl/gsl_sort_vector.h> + + +// Finds the generalized eigenvalues of the Hessian for the state vector z when +// Lambda = c. +void domain_eigen_values(gsl_vector *eigenvalues, unsigned size, unsigned params, gsl_matrix *hess) { + + double eigenvalue; + + gsl_vector *beta; + gsl_vector_complex *alpha; + gsl_matrix *modI; + gsl_eigen_gen_workspace *w; + + alpha = gsl_vector_complex_alloc(size); + beta = gsl_vector_alloc(size); + modI = gsl_matrix_alloc(size, size); + w = gsl_eigen_gen_alloc(size); + + gsl_matrix_set_zero(modI); + for (unsigned i = 0; i < params; i++) gsl_matrix_set(modI, i, i, 1); + + gsl_eigen_gen(hess, modI, alpha, beta, w); + + for (unsigned i = 0; i < size; i++) { + eigenvalue = gsl_vector_complex_get(alpha, i).dat[0] / gsl_vector_get(beta, i); + gsl_vector_set(eigenvalues, i, eigenvalue); + } + + gsl_vector_free(beta); + gsl_vector_complex_free(alpha); + gsl_matrix_free(modI); + gsl_eigen_gen_free(w); +} + + +void domain_eigen_sort(gsl_permutation *eigenorder, unsigned size, unsigned eigen_num, + const gsl_vector *eigenvalues) { + + unsigned ii; + + gsl_vector *abs_eigenvalues; + + abs_eigenvalues = gsl_vector_alloc(size); + + for (unsigned i = 0; i < size; i++) { + gsl_vector_set(abs_eigenvalues, i, fabs(gsl_vector_get(eigenvalues, i))); + } + + gsl_sort_vector_index(eigenorder, abs_eigenvalues); + + gsl_vector_memcpy(abs_eigenvalues, eigenvalues); + + for (unsigned i = eigen_num; i < size; i++) { + ii = gsl_permutation_get(eigenorder, i); + + gsl_vector_set(abs_eigenvalues, ii, INFINITY); + } + + gsl_sort_vector_index(eigenorder, abs_eigenvalues); + + gsl_vector_free(abs_eigenvalues); +} + + +void domain_eigen_state(gsl_vector *eigenstate, const gsl_vector *eigenvalues, + unsigned n, double thres) { + + double eigenvalue; + + gsl_vector_set_zero(eigenstate); + + for (unsigned i = 0; i < 3 * n + 3; i++) { + eigenvalue = gsl_vector_get(eigenvalues, i); + if (eigenvalue > fabs(thres)) gsl_vector_set(eigenstate, i, 1); + if (eigenvalue < -fabs(thres)) gsl_vector_set(eigenstate, i, -1); + } +} + + +void domain_eigen_vector(gsl_vector *eigenvector, unsigned size, unsigned params, unsigned k, gsl_matrix *hess) { + + gsl_vector *beta; + gsl_vector_complex *alpha; + gsl_matrix *modI; + gsl_matrix_complex *evec; + gsl_eigen_genv_workspace *w; + + alpha = gsl_vector_complex_alloc(size); + beta = gsl_vector_alloc(size); + modI = gsl_matrix_alloc(size, size); + evec = gsl_matrix_complex_alloc(size, size); + w = gsl_eigen_genv_alloc(size); + + gsl_matrix_set_zero(modI); + + for (unsigned i = 0; i < params; i++) gsl_matrix_set(modI, i, i, 1); + + gsl_eigen_genv(hess, modI, alpha, beta, evec, w); + + for (unsigned i = 0; i < size; i++) { + gsl_vector_set(eigenvector, i, + gsl_matrix_complex_get(evec, i, k).dat[0]); + } + + gsl_vector_free(beta); + gsl_vector_complex_free(alpha); + gsl_matrix_free(modI); + gsl_matrix_complex_free(evec); + gsl_eigen_genv_free(w); +} + + +void domain_eigen_ortho(gsl_vector *eigenvector, unsigned n, const gsl_vector *z) { + gsl_vector *rotation, *translation_x, *translation_y; + double x, y, prod; + + rotation = gsl_vector_alloc(3 * n + 3); + translation_x = gsl_vector_alloc(3 * n + 3); + translation_y = gsl_vector_alloc(3 * n + 3); + + for (unsigned i = 0; i < n; i++) { + x = gsl_vector_get(z, i); + y = 0; + + gsl_vector_set(translation_x, i, 1.0 / n); + + if (n != 0) { + y = gsl_vector_get(z, n + i - 1); + gsl_vector_set(translation_y, n + i - 1, 1.0 / n); + } + + gsl_vector_set(rotation, i, - y / (gsl_pow_2(x) + gsl_pow_2(y))); + if (n != 0) gsl_vector_set(rotation, n + i - 1, x / (n * (gsl_pow_2(x) + gsl_pow_2(y)))); + } + + gsl_blas_ddot(rotation, eigenvector, &prod); + prod = prod / gsl_blas_dnrm2(rotation); + gsl_vector_memcpy(rotation, eigenvector); + gsl_blas_daxpy(-prod, rotation, eigenvector); + + gsl_blas_ddot(translation_x, eigenvector, &prod); + prod = prod / gsl_blas_dnrm2(translation_x); + gsl_vector_memcpy(translation_x, eigenvector); + gsl_blas_daxpy(-prod, translation_x, eigenvector); + + gsl_blas_ddot(translation_y, eigenvector, &prod); + prod = prod / gsl_blas_dnrm2(translation_y); + gsl_vector_memcpy(translation_y, eigenvector); + gsl_blas_daxpy(-prod, translation_y, eigenvector); + + return; +} + + diff --git a/src/domain_eigen.h b/src/domain_eigen.h new file mode 100644 index 0000000..38366b0 --- /dev/null +++ b/src/domain_eigen.h @@ -0,0 +1,28 @@ +#ifndef DOMAIN_EIGEN_H +#define DOMAIN_EIGEN_H + +#include <gsl/gsl_sf.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_complex.h> +#include <gsl/gsl_complex_math.h> +#include <gsl/gsl_vector.h> +#include <gsl/gsl_permutation.h> +#include <gsl/gsl_permute_vector.h> +#include <gsl/gsl_blas.h> +#include <gsl/gsl_linalg.h> +#include <gsl/gsl_sort_vector.h> + +void domain_eigen_values(gsl_vector *eigenvalues, unsigned size, unsigned params, + gsl_matrix *hess); + +void domain_eigen_sort(gsl_permutation *eigenorder, unsigned size, unsigned eigen_num, + const gsl_vector *eigenvalues); + +void domain_eigen_state(gsl_vector *eigenstate, const gsl_vector *eigenvalues, + unsigned n, double thres); + +void domain_eigen_vector(gsl_vector *eigenvector, unsigned size, unsigned params, unsigned k, gsl_matrix *hess); + +void domain_eigen_ortho(gsl_vector *eigenvector, unsigned n, const gsl_vector *z); + +#endif diff --git a/src/domain_energy.cpp b/src/domain_energy.cpp new file mode 100644 index 0000000..11d6bfe --- /dev/null +++ b/src/domain_energy.cpp @@ -0,0 +1,2103 @@ +/* domain_energy.cpp + * + * Copyright (C) 2013 Jaron Kent-Dobias + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +/* A lightweight and efficient model of two dimensional phase-modulated domains + * using Newton's method to minimize a Lagrangian. + * + * Best viewed in an 80 character terminal with two character hard tabs. + */ + + +#include <stdlib.h> +#include <math.h> + +// GSL includes. +#include <gsl/gsl_sf.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_complex.h> +#include <gsl/gsl_complex_math.h> +#include <gsl/gsl_vector.h> +#include <gsl/gsl_permutation.h> +#include <gsl/gsl_blas.h> + + +void domain_energy_x(gsl_vector *x, unsigned n, const gsl_vector *z) { +// Gets the full set of x coordinates from the state array. + + #pragma omp parallel for + for (unsigned i = 0; i < n; i++) gsl_vector_set(x, i, gsl_vector_get(z, i)); + +} + + +void domain_energy_y(gsl_vector *y, unsigned n, const gsl_vector *z) { +// Gets the full set of y coordinates from the state array. + + gsl_vector_set(y, 0, 0); + + #pragma omp parallel for + for (unsigned i = 1; i < n; i++) { + gsl_vector_set(y, i, gsl_vector_get(z, n + i - 1)); + } + +} + + +void gsl_permutation_over(unsigned n, gsl_permutation *perm, bool right) { +// Shifts a GSL permutation object circularly. If right is true, then the +// permutation is shifted to the right; if false, then it is shifted to the +// left. + + gsl_permutation_swap(perm, 0, n - 1); + + if (right) { + for (unsigned i = 0; i < n - 2; i++) { + gsl_permutation_swap(perm, n - 1 - i, n - 2 - i); + } + } + + else { + for (unsigned i = 0; i < n - 2; i++) { + gsl_permutation_swap(perm, i, i + 1); + } + } +} + + +double domain_energy_area(unsigned n, const gsl_vector *x, const gsl_vector *y) { +// Computes the area of a domain. + + double area, x_i, y_i, x_ii, y_ii; + unsigned ii; + + gsl_permutation *indices_left; + + indices_left = gsl_permutation_alloc(n); + gsl_permutation_init(indices_left); + gsl_permutation_over(n, indices_left, false); + + area = 0; + + for (unsigned i = 0; i < n; i++) { + ii = gsl_permutation_get(indices_left, i); + + x_i = gsl_vector_get(x, i); + y_i = gsl_vector_get(y, i); + x_ii = gsl_vector_get(x, ii); + y_ii = gsl_vector_get(y, ii); + + area += (x_i * y_ii - x_ii * y_i) / 2; + } + + gsl_permutation_free(indices_left); + + return area; +} + + +void domain_energy_rt(gsl_vector *rt, unsigned n, const gsl_vector *x, + double p) { +// Converts x and y coordinates to r_x, r_y, t_x, or t_y, depending on input. + + double x_i, x_ii; + unsigned ii; + + gsl_permutation *indices_left; + + indices_left = gsl_permutation_alloc(n); + gsl_permutation_init(indices_left); + gsl_permutation_over(n, indices_left, false); + + #pragma omp parallel for private(ii, x_i, x_ii) + for (unsigned i = 0; i < n; i++) { + ii = gsl_permutation_get(indices_left, i); + + x_i = gsl_vector_get(x, i); + x_ii = gsl_vector_get(x, ii); + + gsl_vector_set(rt, i, x_ii + p * x_i); + } + + gsl_permutation_free(indices_left); +} + + +void domain_energy_tdots(gsl_matrix *tdots, unsigned n, const gsl_vector *tx, + const gsl_vector *ty) { +// Creates a matrix of dotted tangent vectors. + + gsl_matrix_set_zero(tdots); + + gsl_blas_dger(1, tx, tx, tdots); + gsl_blas_dger(1, ty, ty, tdots); +} + + +void domain_energy_dists(gsl_matrix *dists, unsigned n, const gsl_vector *rx, + const gsl_vector *ry) { +// Creates a matrix of distances between points on the domain boundary.. + + double rx_i, rx_j, ry_i, ry_j; + + #pragma omp parallel for private(rx_i, rx_j, ry_i, ry_j) + for (unsigned i = 0; i < n; i++) { + for (unsigned j = 0; j < n; j++) { + rx_i = gsl_vector_get(rx, i); + rx_j = gsl_vector_get(rx, j); + + ry_i = gsl_vector_get(ry, i); + ry_j = gsl_vector_get(ry, j); + + gsl_matrix_set(dists, i, j, + 2 / sqrt(gsl_pow_2(rx_i - rx_j) + gsl_pow_2(ry_i - ry_j))); + } + } + + #pragma omp parallel for + for (unsigned i = 0; i < n; i++) { + gsl_matrix_set(dists, i, i, 0); + } +} + + +double domain_energy_dconst(unsigned n, double L, double *ld, + const gsl_vector *tx, const gsl_vector *ty) { +// Computes the value of the Lagrangian constraint on the distances. + + double dconst, tx_i, ty_i; + + dconst = 0; + + for (unsigned i = 0; i < n; i++) { + tx_i = gsl_vector_get(tx, i); + ty_i = gsl_vector_get(ty, i); + + dconst += ld[i] * (gsl_pow_2(L / n) + - (gsl_pow_2(tx_i) + gsl_pow_2(ty_i))); + } + + return dconst; +} + + +double domain_energy_dval(unsigned n, double L, double *ld, + const gsl_vector *tx, const gsl_vector *ty) { +// Computes the value of the Lagrangian constraint on the distances. + + double dconst, tx_i, ty_i; + + dconst = 0; + + for (unsigned i = 0; i < n; i++) { + tx_i = gsl_vector_get(tx, i); + ty_i = gsl_vector_get(ty, i); + + dconst += gsl_pow_2(gsl_pow_2(L / n) + - (gsl_pow_2(tx_i) + gsl_pow_2(ty_i))); + } + + return dconst; +} + + +double domain_energy_init(unsigned n, const gsl_vector *z, gsl_vector *rx, + gsl_vector *ry, gsl_vector *tx, gsl_vector *ty, gsl_matrix *tdots, + gsl_matrix *dists) { +// Get useful objects from the state vector. Fills all GSL matrix and vector +// objects, and returns the value of the area. + + double area; + + gsl_vector *x, *y; + + x = gsl_vector_alloc(n); + y = gsl_vector_alloc(n); + + domain_energy_x(x, n, z); + domain_energy_y(y, n, z); + + domain_energy_rt(rx, n, x, 1); + domain_energy_rt(ry, n, y, 1); + domain_energy_rt(tx, n, x, -1); + domain_energy_rt(ty, n, y, -1); + + domain_energy_tdots(tdots, n, tx, ty); + domain_energy_dists(dists, n, rx, ry); + + area = domain_energy_area(n, x, y); + + gsl_vector_free(x); + gsl_vector_free(y); + + return area; +} + + +double domain_energy_energy(unsigned n, + double c, const gsl_vector *rx, const gsl_vector *ry, const gsl_vector *tx, const gsl_vector *ty, gsl_matrix *tdots, gsl_matrix *dists, double L) { +// Computes the Lagrangian. + + double energy_mu, energy, lagrangian, H; + + gsl_vector *v_temp_a, *v_temp_b; + + v_temp_a = gsl_vector_alloc(n); + v_temp_b = gsl_vector_alloc(n); + + gsl_vector_set_all(v_temp_a, 1); + gsl_vector_set_all(v_temp_b, 1); + + gsl_matrix_mul_elements(tdots, dists); + gsl_blas_dtrmv(CblasUpper, CblasNoTrans, CblasNonUnit, tdots, v_temp_a); + gsl_blas_ddot(v_temp_a, v_temp_b, &energy_mu); + + H = log(0.5 * (n - 1)) + M_EULER + 1.0 / (n - 1) + - 1.0 / (12 * gsl_pow_2(0.5 * (n - 1))); + + L = fabs(L); + + energy = c * L - L * log(L) - energy_mu + L * H; + + gsl_vector_free(v_temp_a); + gsl_vector_free(v_temp_b); + + return energy; +} + +double domain_energy_lagrangian(unsigned n, + double c, const gsl_vector *rx, const gsl_vector *ry, const gsl_vector *tx, const gsl_vector *ty, gsl_matrix *tdots, gsl_matrix *dists, double L, double area, double la, double *ld) { +// Computes the Lagrangian. + + double energy, lagrangian; + + energy = domain_energy_energy(n, c, rx, ry, tx, ty, tdots, dists, L); + + L = fabs(L); + + lagrangian = energy - la * (area - M_PI) - domain_energy_dconst(n, L, ld, tx, ty); + + return lagrangian; +} + + +void domain_energy_gradient(gsl_vector *grad, unsigned n, + double c, const gsl_vector *rx, const gsl_vector *ry, const gsl_vector *tx, const gsl_vector *ty, const gsl_matrix *tdots, const gsl_matrix *dists, double L, double area, double la, double *ld) { + +// Computes the gradient of the Lagrangian. + double rx_i, rx_ii, rx_j, tx_i, tx_ii, tx_j, ry_i, ry_ii, + ry_j, ty_i, ty_ii, ty_j, d_ij, d_iij, tdt_ij, tdt_iij, d_ij3, d_iij3; + unsigned ii, jj; + + gsl_vector *v_ones, *v_storage; + gsl_matrix *m_dx, *m_dy; + gsl_permutation *indices_right; + gsl_permutation *indices_left; + + v_ones = gsl_vector_alloc(n); + v_storage = gsl_vector_alloc(n); + + m_dx = gsl_matrix_alloc(n, n); + m_dy = gsl_matrix_alloc(n, n); + + indices_right = gsl_permutation_alloc(n); + indices_left = gsl_permutation_alloc(n); + + gsl_vector_set_zero(grad); + gsl_vector_set_all(v_ones, 1); + gsl_permutation_init(indices_right); + gsl_permutation_over(n, indices_right, true); + gsl_permutation_init(indices_left); + gsl_permutation_over(n, indices_left, false); + + #pragma omp parallel for private(rx_i, rx_j, tx_j, tdt_ij, d_ij, rx_ii, ry_i, ry_j, ty_j, ry_ii, tdt_iij, d_iij, d_ij3, d_iij3, ii) + for (unsigned i = 0; i < n; i++) { + for (unsigned j = 0; j < n; j++) { + ii = gsl_permutation_get(indices_right, i); + + rx_i = gsl_vector_get(rx, i); + rx_ii = gsl_vector_get(rx, ii); + rx_j = gsl_vector_get(rx, j); + tx_j = gsl_vector_get(tx, j); + + ry_i = gsl_vector_get(ry, i); + ry_ii = gsl_vector_get(ry, ii); + ry_j = gsl_vector_get(ry, j); + ty_j = gsl_vector_get(ty, j); + + d_ij = gsl_matrix_get(dists, i, j); + d_iij = gsl_matrix_get(dists, ii, j); + tdt_ij = gsl_matrix_get(tdots, i, j); + tdt_iij = gsl_matrix_get(tdots, ii, j); + + d_ij3 = gsl_pow_3(d_ij); + d_iij3 = gsl_pow_3(d_iij); + + gsl_matrix_set(m_dx, i, j, + - tx_j * d_ij - (rx_i - rx_j) * tdt_ij * d_ij3 / 4 + + tx_j * d_iij - (rx_ii - rx_j) * tdt_iij * d_iij3 / 4); + + gsl_matrix_set(m_dy, i, j, + - ty_j * d_ij - (ry_i - ry_j) * tdt_ij * d_ij3 / 4 + + ty_j * d_iij - (ry_ii - ry_j) * tdt_iij * d_iij3 / 4); + } + } + + #pragma omp parallel for private(ii) + for (unsigned i = 0; i < n; i++) { + ii = gsl_permutation_get(indices_right, i); + gsl_matrix_set(m_dx, i, ii, 0); + gsl_matrix_set(m_dx, i, i, 0); + gsl_matrix_set(m_dy, i, ii, 0); + gsl_matrix_set(m_dy, i, i, 0); + } + + gsl_blas_dgemv(CblasNoTrans, 1, m_dx, v_ones, 0, v_storage); + + #pragma omp parallel for private(ii, tx_i, tx_ii, d_ij) + for (unsigned i = 0; i < n; i++) { + ii = gsl_permutation_get(indices_right, i); + + tx_i = gsl_vector_get(tx, i); + tx_ii = gsl_vector_get(tx, ii); + d_ij = gsl_matrix_get(dists, i, ii); + + gsl_vector_set(grad, i, + - (gsl_vector_get(v_storage, i) + (tx_i - tx_ii) * d_ij)); + } + + gsl_blas_dgemv(CblasNoTrans, 1, m_dy, v_ones, 0, v_storage); + + #pragma omp parallel for private(ii, ty_i, ty_ii, d_ij) + for (unsigned i = 0; i < n; i++) { + ii = gsl_permutation_get(indices_right, i); + + ty_i = gsl_vector_get(ty, i); + ty_ii = gsl_vector_get(ty, ii); + d_ij = gsl_matrix_get(dists, i, ii); + + gsl_vector_set(grad, i + n, + - (gsl_vector_get(v_storage, i) + (ty_i - ty_ii) * d_ij)); + } + + // darea/dx_i or y_i + + #pragma omp parallel for private(ii, jj) + for (unsigned i = 0; i < n; i++) { + ii = gsl_permutation_get(indices_right, i); + jj = gsl_permutation_get(indices_left, i); + + gsl_vector_set(grad, i, gsl_vector_get(grad, i) - + 0.5 * (gsl_vector_get(ty, i) + gsl_vector_get(ty, ii)) * la); + + gsl_vector_set(grad, n + i, gsl_vector_get(grad, n + i) + + 0.5 * (gsl_vector_get(tx, i) + gsl_vector_get(tx, ii)) * la); + } + + #pragma omp parallel for private(ii, jj, tx_i, tx_ii, ty_i, ty_ii) + for (unsigned i = 0; i < n; i++) { + ii = gsl_permutation_get(indices_right, i); + jj = gsl_permutation_get(indices_left, i); + + tx_i = gsl_vector_get(tx, i); + tx_ii = gsl_vector_get(tx, ii); + + ty_i = gsl_vector_get(ty, i); + ty_ii = gsl_vector_get(ty, ii); + + gsl_vector_set(grad, i, gsl_vector_get(grad, i) - + 2 * (ld[i] * tx_i - ld[ii] * tx_ii)); + + gsl_vector_set(grad, i + n, gsl_vector_get(grad, i + n) - + 2 * (ld[i] * ty_i - ld[ii] * ty_ii)); + } + + // The gradient with respect to L. + + L = fabs(L); + double gradLDist = 0; + for (unsigned i = 0; i < n; i++) { + ii = gsl_permutation_get(indices_right, i); + jj = gsl_permutation_get(indices_left, i); + + tx_i = gsl_vector_get(tx, i); + tx_ii = gsl_vector_get(tx, ii); + + ty_i = gsl_vector_get(ty, i); + ty_ii = gsl_vector_get(ty, ii); + + gradLDist += 2 * L / gsl_pow_2(n) * ld[i]; + } + + double H = log(0.5 * (n - 1)) + M_EULER + 1.0 / (n - 1) + - 1.0 / (12 * gsl_pow_2(0.5 * (n - 1))); + + double gradL = GSL_SIGN(L) * (c - (1 + log(L) - H) - gradLDist); + + gsl_vector_set(grad, 2 * n, gradL); + + // The gradients with respect to the undetermined coefficients are simply + // the constraints. + + double gradla = M_PI - area; + + gsl_vector_set(grad, 2 * n + 1, gradla); + + #pragma omp parallel for + for (unsigned i = 0; i < n; i++) { + gsl_vector_set(grad, 2 * n + 2 + i, gsl_pow_2(gsl_vector_get(tx, i)) + + gsl_pow_2(gsl_vector_get(ty, i)) - gsl_pow_2(L / n)); + } + + gsl_vector_free(v_ones); + gsl_vector_free(v_storage); + + gsl_matrix_free(m_dx); + gsl_matrix_free(m_dy); + + gsl_permutation_free(indices_right); + gsl_permutation_free(indices_left); +} + + +void domain_energy_halfHessian(gsl_matrix *hess, unsigned n, + double c, const gsl_vector *rx, const gsl_vector *ry, const gsl_vector *tx, const gsl_vector *ty, const gsl_matrix *tdots, const gsl_matrix *dists, double L, double la, double *ld) { +/* Computes the Hessian of the Lagrangian without the center of mass + * constraints and fixed point. + */ + + gsl_matrix *m_dxidxj, *m_dyidyj, *m_dxidxi, *m_dyidyi, *m_dxidxii, + *m_dyidyii, *m_dxidyj, *m_dxidyi, *m_dxiidyi, *m_dxidyii; + gsl_vector *v_ones, *v_storage; + gsl_permutation *indicesRight, *indicesLeft, *indices2Right; + + unsigned ii, jj, iii; + double rx_i, rx_j, tx_i, tx_j, tdt_ij, d_ij, rx_ii, rx_jj, tx_ii, tx_jj, + ry_i, ry_j, ty_i, ty_j, ry_ii, ry_jj, ty_ii, ty_jj, tdt_iij, d_iij, + tdt_ijj, d_ijj, tdt_iijj, d_iijj, rx_iii, tx_iii, ry_iii, ty_iii, + d_ij3, d_iij3, d_ijj3, d_iijj3, d_ij5, d_iij5, d_ijj5, d_iijj5; + + m_dxidxj = gsl_matrix_alloc(n, n); + m_dyidyj = gsl_matrix_alloc(n, n); + m_dxidxi = gsl_matrix_alloc(n, n); + m_dyidyi = gsl_matrix_alloc(n, n); + m_dxidxii = gsl_matrix_alloc(n, n); + m_dyidyii = gsl_matrix_alloc(n, n); + m_dxidyj = gsl_matrix_alloc(n, n); + m_dxidyi = gsl_matrix_alloc(n, n); + m_dxiidyi = gsl_matrix_alloc(n, n); + m_dxidyii = gsl_matrix_alloc(n, n); + + v_ones = gsl_vector_alloc(n); + v_storage = gsl_vector_alloc(n); + + indicesRight = gsl_permutation_alloc(n); + indicesLeft = gsl_permutation_alloc(n); + indices2Right = gsl_permutation_alloc(n); + + gsl_matrix_set_zero(hess); + gsl_vector_set_all(v_ones, 1); + + gsl_permutation_init(indicesRight); + gsl_permutation_init(indicesLeft); + gsl_permutation_over(n, indicesRight, true); + gsl_permutation_over(n, indicesLeft, false); + gsl_permutation_memcpy(indices2Right, indicesRight); + gsl_permutation_over(n, indices2Right, true); + + #pragma omp parallel for private(rx_i, rx_j, tx_i, tx_j, tdt_ij, d_ij, rx_ii, rx_jj, tx_ii, tx_jj, ry_i, ry_j, ty_i, ty_j, ry_ii, ry_jj, ty_ii, ty_jj, tdt_iij, d_iij, tdt_ijj, d_ijj, tdt_iijj, d_iijj, rx_iii, tx_iii, ry_iii, ty_iii, d_ij3, d_iij3, d_ijj3, d_iijj3, d_ij5, d_iij5, d_ijj5, d_iijj5, ii, jj) + for (unsigned i = 0; i < n; i++) { + ii = gsl_permutation_get(indicesRight, i); + for (unsigned j = 0; j < n; j++) { + jj = gsl_permutation_get(indicesRight, j); + + rx_i = gsl_vector_get(rx, i); + rx_j = gsl_vector_get(rx, j); + tx_i = gsl_vector_get(tx, i); + tx_j = gsl_vector_get(tx, j); + + ry_i = gsl_vector_get(ry, i); + ry_j = gsl_vector_get(ry, j); + ty_i = gsl_vector_get(ty, i); + ty_j = gsl_vector_get(ty, j); + + rx_ii = gsl_vector_get(rx, ii); + rx_jj = gsl_vector_get(rx, jj); + tx_ii = gsl_vector_get(tx, ii); + tx_jj = gsl_vector_get(tx, jj); + + ry_ii = gsl_vector_get(ry, ii); + ry_jj = gsl_vector_get(ry, jj); + ty_ii = gsl_vector_get(ty, ii); + ty_jj = gsl_vector_get(ty, jj); + + d_ij = gsl_matrix_get(dists, i, j); + tdt_ij = gsl_matrix_get(tdots, i, j); + + d_iij = gsl_matrix_get(dists, ii, j); + tdt_iij = gsl_matrix_get(tdots, ii, j); + + d_ijj = gsl_matrix_get(dists, i, jj); + tdt_ijj = gsl_matrix_get(tdots, i, jj); + + d_iijj = gsl_matrix_get(dists, ii, jj); + tdt_iijj = gsl_matrix_get(tdots, ii, jj); + + d_ij3 = gsl_pow_3(d_ij); + d_iij3 = gsl_pow_3(d_iij); + d_ijj3 = gsl_pow_3(d_ijj); + d_iijj3 = gsl_pow_3(d_iijj); + d_ij5 = gsl_pow_5(d_ij); + d_iij5 = gsl_pow_5(d_iij); + d_ijj5 = gsl_pow_5(d_ijj); + d_iijj5 = gsl_pow_5(d_iijj); + + // d^2E / dx_i dx_j for i != j, j - 1. + gsl_matrix_set(m_dxidxj, i, j, + (rx_i - rx_j) * (tx_i - tx_j) * d_ij3 / 4 + d_ij + - 3 * gsl_pow_2(rx_i - rx_j) * tdt_ij * d_ij5 / 16 + + tdt_ij * d_ij3 / 4 + + + (rx_ii - rx_j) * (tx_ii + tx_j) * d_iij3 / 4 - d_iij + - 3 * gsl_pow_2(rx_ii - rx_j) * tdt_iij * d_iij5 / 16 + + tdt_iij * d_iij3 / 4 + + - (rx_i - rx_jj) * (tx_i + tx_jj) * d_ijj3 / 4 - d_ijj + - 3 * gsl_pow_2(rx_i - rx_jj) * tdt_ijj * d_ijj5 / 16 + + tdt_ijj * d_ijj3 / 4 + + - (rx_ii - rx_jj) * (tx_ii - tx_jj) * d_iijj3 / 4 + d_iijj + - 3 * gsl_pow_2(rx_ii - rx_jj) * tdt_iijj * d_iijj5 / 16 + + tdt_iijj * d_iijj3 / 4 + ); + + // d^2E / dy_i dy_j for i != j, j - 1. + gsl_matrix_set(m_dyidyj, i, j, + (ry_i - ry_j) * (ty_i - ty_j) * d_ij3 / 4 + d_ij + - 3 * gsl_pow_2(ry_i - ry_j) * tdt_ij * d_ij5 / 16 + + tdt_ij * d_ij3 / 4 + + + (ry_ii - ry_j) * (ty_ii + ty_j) * d_iij3 / 4 - d_iij + - 3 * gsl_pow_2(ry_ii - ry_j) * tdt_iij * d_iij5 / 16 + + tdt_iij * d_iij3 / 4 + + - (ry_i - ry_jj) * (ty_i + ty_jj) * d_ijj3 / 4 - d_ijj + - 3 * gsl_pow_2(ry_i - ry_jj) * tdt_ijj * d_ijj5 / 16 + + tdt_ijj * d_ijj3 / 4 + + - (ry_ii - ry_jj) * (ty_ii - ty_jj) * d_iijj3 / 4 + d_iijj + - 3 * gsl_pow_2(ry_ii - ry_jj) * tdt_iijj * d_iijj5 / 16 + + tdt_iijj * d_iijj3 / 4 + ); + + // d^2E / dx_i^2 + gsl_matrix_set(m_dxidxi, i, j, + (rx_i - rx_j) * tx_j * d_ij3 / 2 + + 3 * gsl_pow_2(rx_i - rx_j) * tdt_ij * d_ij5 / 16 + - tdt_ij * d_ij3 / 4 + + - (rx_ii - rx_j) * tx_j * d_iij3 / 2 + + 3 * gsl_pow_2(rx_ii - rx_j) * tdt_iij * d_iij5 / 16 + - tdt_iij * d_iij3 / 4 + ); + + // d^2E / dy_i^2 + gsl_matrix_set(m_dyidyi, i, j, + (ry_i - ry_j) * ty_j * d_ij3 / 2 + + 3 * gsl_pow_2(ry_i - ry_j) * tdt_ij * d_ij5 / 16 + - tdt_ij * d_ij3 / 4 + + - (ry_ii - ry_j) * ty_j * d_iij3 / 2 + + 3 * gsl_pow_2(ry_ii - ry_j) * tdt_iij * d_iij5 / 16 + - tdt_iij * d_iij3 / 4 + ); + + // d^2E / dx_i dx_(i-1) + gsl_matrix_set(m_dxidxii, i, j, + 3 * gsl_pow_2(rx_ii - rx_j) * tdt_iij * d_iij5 / 16 + - tdt_iij * d_iij3 / 4 + ); + + // d^2E / dy_i dy_(i-1) + gsl_matrix_set(m_dyidyii, i, j, + 3 * gsl_pow_2(ry_ii - ry_j) * tdt_iij * d_iij5 / 16 + - tdt_iij * d_iij3 / 4 + ); + + gsl_matrix_set(m_dxidyj, i, j, + (rx_i - rx_j) * ty_i * d_ij3 / 4 + - (ry_i - ry_j) * tx_j * d_ij3 / 4 + - 3 * (rx_i - rx_j) * (ry_i - ry_j) * tdt_ij * d_ij5 / 16 + + + (rx_ii - rx_j) * ty_ii * d_iij3 / 4 + + (ry_ii - ry_j) * tx_j * d_iij3 / 4 + - 3 * (rx_ii - rx_j) * (ry_ii - ry_j) * tdt_iij * d_iij5 / 16 + + - (rx_i - rx_jj) * ty_i * d_ijj3 / 4 + - (ry_i - ry_jj) * tx_jj * d_ijj3 / 4 + - 3 * (rx_i - rx_jj) * (ry_i - ry_jj) * tdt_ijj * d_ijj5 / 16 + + - (rx_ii - rx_jj) * ty_ii * d_iijj3 / 4 + + (ry_ii - ry_jj) * tx_jj * d_iijj3 / 4 + - 3 * (rx_ii - rx_jj) * (ry_ii - ry_jj) * tdt_iijj * d_iijj5 / 16 + ); + + gsl_matrix_set(m_dxidyi, i, j, + (ry_i - ry_j) * tx_j * d_ij3 / 4 + + (rx_i - rx_j) * ty_j * d_ij3 / 4 + + 3 * (rx_i - rx_j) * (ry_i - ry_j) * tdt_ij * d_ij5 / 16 + + - (ry_ii - ry_j) * tx_j * d_iij3 / 4 + - (rx_ii - rx_j) * ty_j * d_iij3 / 4 + + 3 * (rx_ii - rx_j) * (ry_ii - ry_j) * tdt_iij * d_iij5 / 16 + ); + + gsl_matrix_set(m_dxiidyi, i, j, + (ry_ii - ry_j) * tx_j * d_iij3 / 4 + - (rx_ii - rx_j) * ty_j * d_iij3 / 4 + + 3 * (rx_ii - rx_j) * (ry_ii - ry_j) * tdt_iij * d_iij5 / 16 + ); + + gsl_matrix_set(m_dxidyii, i, j, + - (ry_i - ry_j) * tx_j * d_ij3 / 4 + + (rx_i - rx_j) * ty_j * d_ij3 / 4 + + 3 * (rx_i - rx_j) * (ry_i - ry_j) * tdt_ij * d_ij5 / 16 + ); + + } + } + + // Setting terms of d^2E / dy_i dy_j and d^2E / dx_i dx_j in the Hessian. + #pragma omp parallel for + for (unsigned i = 2; i < n; i++) { + for (unsigned j = 0; j < i - 1; j++) { + + gsl_matrix_set(hess, i, j, - gsl_matrix_get(m_dxidxj, i, j)); + + gsl_matrix_set(hess, n + i, n + j, - gsl_matrix_get(m_dyidyj, i, j)); + } + } + + // Zeroing out terms which aren't supposed to appear in the sums. + #pragma omp parallel for private(ii, iii, jj) + for (unsigned i = 0; i < n; i++) { + ii = gsl_permutation_get(indicesRight, i); + iii = gsl_permutation_get(indices2Right, i); + jj = gsl_permutation_get(indicesLeft, i); + gsl_matrix_set(m_dxidxi, i, i, 0); + gsl_matrix_set(m_dxidxi, i, ii, 0); + gsl_matrix_set(m_dyidyi, i, i, 0); + gsl_matrix_set(m_dyidyi, i, ii, 0); + gsl_matrix_set(m_dxidxii, i, i, 0); + gsl_matrix_set(m_dxidxii, i, ii, 0); + gsl_matrix_set(m_dxidxii, i, iii, 0); + gsl_matrix_set(m_dyidyii, i, i, 0); + gsl_matrix_set(m_dyidyii, i, ii, 0); + gsl_matrix_set(m_dyidyii, i, iii, 0); + gsl_matrix_set(m_dxidyi, i, i, 0); + gsl_matrix_set(m_dxidyi, i, ii, 0); + gsl_matrix_set(m_dxiidyi, i, i, 0); + gsl_matrix_set(m_dxiidyi, i, ii, 0); + gsl_matrix_set(m_dxiidyi, i, iii, 0); + gsl_matrix_set(m_dxidyii, i, ii, 0); + gsl_matrix_set(m_dxidyii, i, jj, 0); + gsl_matrix_set(m_dxidyii, i, i, 0); + } + + gsl_blas_dgemv(CblasNoTrans, 1, m_dxidxi, v_ones, 0, v_storage); + + // Setting terms of d^2E / dx_i^2 in the Hessian. + #pragma omp parallel for private(ii, d_iij) + for (unsigned i = 0; i < n; i++) { + ii = gsl_permutation_get(indicesRight, i); + + d_iij = gsl_matrix_get(dists, i, ii); + + gsl_matrix_set(hess, i, i, + - (gsl_vector_get(v_storage, i) - 2 * d_iij) + + 2 * (ld[i] + ld[ii]) + ); + } + + gsl_blas_dgemv(CblasNoTrans, 1, m_dyidyi, v_ones, 0, v_storage); + + // Setting terms of d^2E / dy_i^2 in the Hessian. + #pragma omp parallel for private(ii, d_iij) + for (unsigned i = 0; i < n; i++) { + ii = gsl_permutation_get(indicesRight, i); + + d_iij = gsl_matrix_get(dists, i, ii); + + gsl_matrix_set(hess, n + i, n + i, + - (gsl_vector_get(v_storage, i) - 2 * d_iij) + + 2 * (ld[i] + ld[ii]) + ); + } + + gsl_blas_dgemv(CblasNoTrans, 1, m_dxidxii, v_ones, 0, v_storage); + + #pragma omp parallel for private(ii, rx_i, rx_ii, tx_i, tx_ii, d_ij) + for (unsigned i = 0; i < n; i++) { + ii = gsl_permutation_get(indicesRight, i); + + rx_i = gsl_vector_get(rx, i); + rx_ii = gsl_vector_get(rx, ii); + tx_i = gsl_vector_get(tx, i); + tx_ii = gsl_vector_get(tx, ii); + d_ij = gsl_matrix_get(dists, i, ii); + gsl_vector_set(v_storage, i, gsl_vector_get(v_storage, i) + - (rx_ii - rx_i) * (tx_i - tx_ii) * gsl_pow_3(d_ij) / 4 + + d_ij); + } + + #pragma omp parallel for private(ii, iii, rx_i, rx_ii, rx_iii, tx_i, tx_ii, tx_iii, d_ij, d_iij, d_iijj, tdt_ij) + for (unsigned i = 0; i < n; i++) { + iii = gsl_permutation_get(indices2Right, i); + ii = gsl_permutation_get(indicesRight, i); + + rx_i = gsl_vector_get(rx, i); + rx_ii = gsl_vector_get(rx, ii); + rx_iii = gsl_vector_get(rx, iii); + tx_i = gsl_vector_get(tx, i); + tx_ii = gsl_vector_get(tx, ii); + tx_iii = gsl_vector_get(tx, iii); + d_ij = gsl_matrix_get(dists, i, ii); + d_iij = gsl_matrix_get(dists, i, iii); + d_iijj = gsl_matrix_get(dists, ii, iii); + tdt_ij = gsl_matrix_get(tdots, i, iii); + + gsl_vector_set(v_storage, i, gsl_vector_get(v_storage, i) + - 3 * gsl_pow_2(rx_iii - rx_i) * tdt_ij * gsl_pow_5(d_iij) / 16 + + (tx_iii + tx_i) * (rx_iii - rx_i) * gsl_pow_3(d_iij) / 4 + + tdt_ij * gsl_pow_3(d_iij) / 4 - d_iij + - (tx_ii - tx_iii) * (rx_ii - rx_iii) * gsl_pow_3(d_iijj) / 4 + d_iijj); + } + + gsl_matrix_set(hess, n - 1, 0, + - gsl_vector_get(v_storage, 0) - 2 * ld[n - 1]); + + #pragma omp parallel for + for (unsigned i = 1; i < n; i++) { + gsl_matrix_set(hess, i, i - 1, + - gsl_vector_get(v_storage, i) - 2 * ld[i - 1] + ); + } + + + gsl_blas_dgemv(CblasNoTrans, 1, m_dyidyii, v_ones, 0, v_storage); + + #pragma omp parallel for private(ii, ry_i, ry_ii, ty_i, ty_ii, d_ij) + for (unsigned i = 0; i < n; i++) { + ii = gsl_permutation_get(indicesRight, i); + + ry_i = gsl_vector_get(ry, i); + ry_ii = gsl_vector_get(ry, ii); + ty_i = gsl_vector_get(ty, i); + ty_ii = gsl_vector_get(ty, ii); + d_ij = gsl_matrix_get(dists, i, ii); + gsl_vector_set(v_storage, i, gsl_vector_get(v_storage, i) + - (ry_ii - ry_i) * (ty_i - ty_ii) * gsl_pow_3(d_ij) / 4 + + d_ij); + } + + #pragma omp parallel for private(ii, iii, ry_i, ry_ii, ry_iii, ty_i, ty_ii, ty_iii, d_ij, d_iij, d_iijj, tdt_ij) + for (unsigned i = 0; i < n; i++) { + iii = gsl_permutation_get(indices2Right, i); + ii = gsl_permutation_get(indicesRight, i); + + ry_i = gsl_vector_get(ry, i); + ry_ii = gsl_vector_get(ry, ii); + ry_iii = gsl_vector_get(ry, iii); + ty_i = gsl_vector_get(ty, i); + ty_ii = gsl_vector_get(ty, ii); + ty_iii = gsl_vector_get(ty, iii); + d_ij = gsl_matrix_get(dists, i, ii); + d_iij = gsl_matrix_get(dists, i, iii); + d_iijj = gsl_matrix_get(dists, ii, iii); + tdt_ij = gsl_matrix_get(tdots, i, iii); + + gsl_vector_set(v_storage, i, gsl_vector_get(v_storage, i) + - 3 * gsl_pow_2(ry_iii - ry_i) * tdt_ij * gsl_pow_5(d_iij) / 16 + + (ty_iii + ty_i) * (ry_iii - ry_i) * gsl_pow_3(d_iij) / 4 + + tdt_ij * gsl_pow_3(d_iij) / 4 - d_iij + - (ty_ii - ty_iii) * (ry_ii - ry_iii) * gsl_pow_3(d_iijj) / 4 + d_iijj); + } + + #pragma omp parallel for + for (unsigned i = 1; i < n; i++) { + gsl_matrix_set(hess, n + i, n + i - 1, + - gsl_vector_get(v_storage, i) - 2 * ld[i - 1] + ); + } + + gsl_matrix_set(hess, 2 * n -1, n, + - gsl_vector_get(v_storage, 0) - 2 * ld[n-1]); + + + // dxdy boring style + #pragma omp parallel for + for (unsigned j = 2; j < n; j++) { + for (unsigned i = 0; i < j - 1; i++) { + gsl_matrix_set(hess, n + j, i, + - gsl_matrix_get(m_dxidyj, i, j)); + } + } + + #pragma omp parallel for + for (unsigned j = 0; j < n - 2; j++) { + for (unsigned i = j + 2; i < n; i++) { + gsl_matrix_set(hess, n + j, i, + - gsl_matrix_get(m_dxidyj, i, j)); + } + } + + // d^2E / dx_i dy_i + + gsl_blas_dgemv(CblasNoTrans, 1, m_dxidyi, v_ones, 0, v_storage); + + #pragma omp parallel for + for (unsigned i = 0; i < n; i++) { + gsl_matrix_set(hess, n + i, i, + - gsl_vector_get(v_storage, i) + ); + } + + // d^2 E / dx_ii dy_i + + gsl_blas_dgemv(CblasNoTrans, 1, m_dxiidyi, v_ones, 0, v_storage); + + #pragma omp parallel for private(ii, iii, rx_i, rx_ii, rx_iii, tx_i, tx_ii, tx_iii, ry_i, ry_ii, ry_iii, ty_i, ty_ii, ty_iii, d_ij, d_iij, d_iijj, tdt_iij) + for (unsigned i = 0; i < n; i++) { + ii = gsl_permutation_get(indicesRight, i); + iii = gsl_permutation_get(indices2Right, i); + + rx_i = gsl_vector_get(rx, i); + rx_ii = gsl_vector_get(rx, ii); + rx_iii = gsl_vector_get(rx, iii); + tx_i = gsl_vector_get(tx, i); + tx_ii = gsl_vector_get(tx, ii); + tx_iii = gsl_vector_get(tx, iii); + ry_i = gsl_vector_get(ry, i); + ry_ii = gsl_vector_get(ry, ii); + ry_iii = gsl_vector_get(ry, iii); + ty_i = gsl_vector_get(ty, i); + ty_ii = gsl_vector_get(ty, ii); + ty_iii = gsl_vector_get(ty, iii); + + d_ij = gsl_matrix_get(dists, i, ii); + d_iij = gsl_matrix_get(dists, i, iii); + d_iijj = gsl_matrix_get(dists, ii, iii); + tdt_iij = gsl_matrix_get(tdots, i, iii); + + gsl_vector_set(v_storage, i, gsl_vector_get(v_storage, i) + - (rx_ii - rx_i) * (ty_i - ty_ii) * gsl_pow_3(d_ij) / 4 + - 3 * (rx_iii - rx_i) * (ry_iii - ry_i) * tdt_iij * gsl_pow_5(d_iij) / 16 + + (rx_iii - rx_i) * ty_iii * gsl_pow_3(d_iij) / 4 + + (ry_iii - ry_i) * tx_i * gsl_pow_3(d_iij) / 4 + - (tx_ii - tx_iii) * (ry_ii - ry_iii) * gsl_pow_3(d_iijj) / 4 + ); + } + + #pragma omp parallel for + for (unsigned i = 0; i < n - 1; i++) { + gsl_matrix_set(hess, n + i + 1, i, + - gsl_vector_get(v_storage, i + 1) - la / 2 + ); + } + + gsl_matrix_set(hess, n, n - 1, + -gsl_vector_get(v_storage, 0) - la / 2); + + // Upper off-diagonal of dxdy submatrix. + + gsl_blas_dgemv(CblasNoTrans, 1, m_dxidyii, v_ones, 0, v_storage); + + #pragma omp parallel for private(ii, jj, rx_i, rx_ii, rx_jj, tx_i, tx_ii, tx_jj, ry_i, ry_ii, ry_jj, ty_i, ty_ii, ty_jj, d_ij, d_ijj, d_iijj, tdt_iijj) + for (unsigned i = 0; i < n; i++) { + ii = gsl_permutation_get(indicesRight, i); + jj = gsl_permutation_get(indicesLeft, i); + + rx_i = gsl_vector_get(rx, i); + rx_ii = gsl_vector_get(rx, ii); + rx_jj = gsl_vector_get(rx, jj); + tx_i = gsl_vector_get(tx, i); + tx_ii = gsl_vector_get(tx, ii); + tx_jj = gsl_vector_get(tx, jj); + ry_i = gsl_vector_get(ry, i); + ry_ii = gsl_vector_get(ry, ii); + ry_jj = gsl_vector_get(ry, jj); + ty_i = gsl_vector_get(ty, i); + ty_ii = gsl_vector_get(ty, ii); + ty_jj = gsl_vector_get(ty, jj); + + d_ij = gsl_matrix_get(dists, i, ii); + d_ijj = gsl_matrix_get(dists, i, jj); + d_iijj = gsl_matrix_get(dists, ii, jj); + tdt_iijj = gsl_matrix_get(tdots, ii, jj); + + gsl_vector_set(v_storage, i, gsl_vector_get(v_storage, i) + - (rx_i - rx_ii) * (ty_i - ty_ii) * gsl_pow_3(d_ij) / 4 + - (tx_jj - tx_i) * (ry_i - ry_jj) * gsl_pow_3(d_ijj) / 4 + + tx_ii * (ry_ii - ry_jj) * gsl_pow_3(d_iijj) / 4 + - (rx_jj - rx_ii) * ty_jj * gsl_pow_3(d_iijj) / 4 + + 3 * (rx_jj - rx_ii) * (ry_ii - ry_jj) * tdt_iijj * gsl_pow_5(d_iijj) / 16 + ); + } + + gsl_matrix_set(hess, 2 * n - 1, 0, + - gsl_vector_get(v_storage, n - 1) + la / 2 + ); + + + #pragma omp parallel for + for (unsigned i = 0; i < n-1; i++) { + gsl_matrix_set(hess, n + i, i + 1, + - gsl_vector_get(v_storage, i) + la / 2 + ); + } + + + // dLdL + + double gradLDist = 0; + for (unsigned i = 0; i < n; i++) { + gradLDist += 2 * ld[i] / gsl_pow_2(n); + } + + L = fabs(L); + + double gradL = - 1 / L - gradLDist; + + gsl_matrix_set(hess, 2 * n, 2 * n, gradL); + + // dLdlambdad + + #pragma omp parallel for + for (unsigned i = 0; i < n; i++) { + gsl_matrix_set(hess, i + 2 * n + 2, 2 * n, - 2 * L / gsl_pow_2(n)); + } + + #pragma omp parallel for private(ii) + for (unsigned i = 0; i < n; i++) { + ii = gsl_permutation_get(indicesRight, i); + + gsl_matrix_set(hess, 2 * n + 1, i, - + 0.5 * (gsl_vector_get(ty, i) + gsl_vector_get(ty, ii))); + + gsl_matrix_set(hess, 2 * n + 1, n + i, + 0.5 * (gsl_vector_get(tx, i) + gsl_vector_get(tx, ii))); + } + + #pragma omp parallel for private(ii) + for (unsigned i = 0; i < n; i++) { + ii = gsl_permutation_get(indicesRight, i); + + gsl_matrix_set(hess, 2 * n + 2 + i, i, - + 2 * (gsl_vector_get(tx, i))); + + gsl_matrix_set(hess, 2 * n + 2 + i, n + i, - + 2 * (gsl_vector_get(ty, i))); + } + + #pragma omp parallel for private(ii) + for (unsigned i = 1; i < n; i++) { + ii = gsl_permutation_get(indicesRight, i); + + gsl_matrix_set(hess, 2 * n + i + 1, i, + 2 * gsl_vector_get(tx, ii)); + + gsl_matrix_set(hess, 2 * n + i + 1, n + i, + 2 * gsl_vector_get(ty, ii)); + } + + gsl_matrix_set(hess, 3 * n + 1, 0, + 2 * gsl_vector_get(tx, n-1)); + gsl_matrix_set(hess, 3 * n + 1, n, + 2 * gsl_vector_get(ty, n-1)); + + + gsl_vector_free(v_ones); + gsl_vector_free(v_storage); + + gsl_matrix_free(m_dxidxj); + gsl_matrix_free(m_dyidyj); + gsl_matrix_free(m_dxidxi); + gsl_matrix_free(m_dyidyi); + gsl_matrix_free(m_dxidxii); + gsl_matrix_free(m_dyidyii); + gsl_matrix_free(m_dxidyj); + gsl_matrix_free(m_dxidyi); + gsl_matrix_free(m_dxiidyi); + gsl_matrix_free(m_dxidyii); + + gsl_permutation_free(indicesLeft); + gsl_permutation_free(indicesRight); + gsl_permutation_free(indices2Right); +} + + +double domain_energy_nakedEnergy(unsigned n, const gsl_vector *z, double c) { + double lagrangian; + + gsl_vector *x, *y; + + x = gsl_vector_alloc(n); + y = gsl_vector_alloc(n); + + #pragma omp parallel for + for (unsigned i = 0; i < n; i++) gsl_vector_set(x, i, gsl_vector_get(z, i)); + #pragma omp parallel for + for (unsigned i = 0; i < n; i++) gsl_vector_set(y, i, gsl_vector_get(z, i + n)); + double L = gsl_vector_get(z, 2 * n); + + unsigned ii; + + gsl_vector *rx, *ry, *tx, *ty; + gsl_matrix *tdots, *dists; + + rx = gsl_vector_alloc(n); + ry = gsl_vector_alloc(n); + tx = gsl_vector_alloc(n); + ty = gsl_vector_alloc(n); + + tdots = gsl_matrix_alloc(n, n); + dists = gsl_matrix_alloc(n, n); + + domain_energy_rt(rx, n, x, 1); + domain_energy_rt(ry, n, y, 1); + domain_energy_rt(tx, n, x, -1); + domain_energy_rt(ty, n, y, -1); + + domain_energy_tdots(tdots, n, tx, ty); + domain_energy_dists(dists, n, rx, ry); + + lagrangian = domain_energy_energy(n, c, rx, ry, tx, ty, tdots, dists, L); + + gsl_vector_free(x); + gsl_vector_free(y); + gsl_vector_free(rx); + gsl_vector_free(ry); + gsl_vector_free(tx); + gsl_vector_free(ty); + + gsl_matrix_free(tdots); + gsl_matrix_free(dists); + + return lagrangian; +}; + + +double domain_energy_nakedLagrangian(unsigned n, const gsl_vector *z, double c) { + double lagrangian; + + gsl_vector *x, *y; + + x = gsl_vector_alloc(n); + y = gsl_vector_alloc(n); + + #pragma omp parallel for + for (unsigned i = 0; i < n; i++) gsl_vector_set(x, i, gsl_vector_get(z, i)); + #pragma omp parallel for + for (unsigned i = 0; i < n; i++) gsl_vector_set(y, i, gsl_vector_get(z, i + n)); + double L = gsl_vector_get(z, 2 * n); + double la = gsl_vector_get(z, 2 * n + 1); + double ld[n]; + #pragma omp parallel for + for (unsigned i = 0; i < n; i++) { + ld[i] = gsl_vector_get(z, 2 * n + 2 + i); + } + + unsigned ii; + + gsl_vector *rx, *ry, *tx, *ty; + gsl_matrix *tdots, *dists; + + rx = gsl_vector_alloc(n); + ry = gsl_vector_alloc(n); + tx = gsl_vector_alloc(n); + ty = gsl_vector_alloc(n); + + tdots = gsl_matrix_alloc(n, n); + dists = gsl_matrix_alloc(n, n); + + domain_energy_rt(rx, n, x, 1); + domain_energy_rt(ry, n, y, 1); + domain_energy_rt(tx, n, x, -1); + domain_energy_rt(ty, n, y, -1); + + domain_energy_tdots(tdots, n, tx, ty); + domain_energy_dists(dists, n, rx, ry); + + double area = domain_energy_area(n, x, y); + + lagrangian = domain_energy_lagrangian(n, c, rx, ry, tx, ty, tdots, dists, L, area, la, ld); + + gsl_vector_free(x); + gsl_vector_free(y); + gsl_vector_free(rx); + gsl_vector_free(ry); + gsl_vector_free(tx); + gsl_vector_free(ty); + + gsl_matrix_free(tdots); + gsl_matrix_free(dists); + + return lagrangian; +}; + + +void domain_energy_nakedGradient(gsl_vector *grad, unsigned n, const gsl_vector *z, double c) { + + gsl_permutation *indices_right; + indices_right = gsl_permutation_alloc(n); + + gsl_permutation_init(indices_right); + gsl_permutation_over(n, indices_right, true); + + gsl_vector *x, *y; + + x = gsl_vector_alloc(n); + y = gsl_vector_alloc(n); + + // Setting pointers to give the elements of z more convenient names. + #pragma omp parallel for + for (unsigned i = 0; i < n; i++) gsl_vector_set(x, i, gsl_vector_get(z, i)); + #pragma omp parallel for + for (unsigned i = 0; i < n; i++) gsl_vector_set(y, i, gsl_vector_get(z, i + n)); + double L = gsl_vector_get(z, 2 * n); + double la = gsl_vector_get(z, 2 * n + 1); + double ld[n]; + #pragma omp parallel for + for (unsigned i = 0; i < n; i++) { + ld[i] = gsl_vector_get(z, 2 * n + 2 + i); + } + + unsigned ii; + + gsl_vector *rx, *ry, *tx, *ty; + gsl_matrix *tdots, *dists; + + rx = gsl_vector_alloc(n); + ry = gsl_vector_alloc(n); + tx = gsl_vector_alloc(n); + ty = gsl_vector_alloc(n); + + tdots = gsl_matrix_alloc(n, n); + dists = gsl_matrix_alloc(n, n); + + domain_energy_rt(rx, n, x, 1); + domain_energy_rt(ry, n, y, 1); + domain_energy_rt(tx, n, x, -1); + domain_energy_rt(ty, n, y, -1); + + domain_energy_tdots(tdots, n, tx, ty); + domain_energy_dists(dists, n, rx, ry); + + double area = domain_energy_area(n, x, y); + + domain_energy_gradient(grad, n, c, rx, ry, tx, ty, tdots, dists, L, area, la, ld); + + gsl_vector_free(x); + gsl_vector_free(y); + gsl_vector_free(rx); + gsl_vector_free(ry); + gsl_vector_free(tx); + gsl_vector_free(ty); + gsl_matrix_free(tdots); + gsl_matrix_free(dists); + gsl_permutation_free(indices_right); + +} + + +void domain_energy_nakedHalfHessian(gsl_matrix *hess, unsigned n, const gsl_vector *z, double c) { + + gsl_permutation *indices_right; + indices_right = gsl_permutation_alloc(n); + + gsl_permutation_init(indices_right); + gsl_permutation_over(n, indices_right, true); + + gsl_vector *x, *y; + + x = gsl_vector_alloc(n); + y = gsl_vector_alloc(n); + + // Setting pointers to give the elements of z more convenient names. + #pragma omp parallel for + for (unsigned i = 0; i < n; i++) gsl_vector_set(x, i, gsl_vector_get(z, i)); + #pragma omp parallel for + for (unsigned i = 0; i < n; i++) gsl_vector_set(y, i, gsl_vector_get(z, i + n)); + double L = gsl_vector_get(z, 2 * n); + double la = gsl_vector_get(z, 2 * n + 1); + double ld[n]; + #pragma omp parallel for + for (unsigned i = 0; i < n; i++) { + ld[i] = gsl_vector_get(z, 2 * n + 2 + i); + } + + unsigned ii; + + gsl_vector *rx, *ry, *tx, *ty; + gsl_matrix *tdots, *dists; + + rx = gsl_vector_alloc(n); + ry = gsl_vector_alloc(n); + tx = gsl_vector_alloc(n); + ty = gsl_vector_alloc(n); + + tdots = gsl_matrix_alloc(n, n); + dists = gsl_matrix_alloc(n, n); + + domain_energy_rt(rx, n, x, 1); + domain_energy_rt(ry, n, y, 1); + domain_energy_rt(tx, n, x, -1); + domain_energy_rt(ty, n, y, -1); + + domain_energy_tdots(tdots, n, tx, ty); + domain_energy_dists(dists, n, rx, ry); + + double area = domain_energy_area(n, x, y); + + domain_energy_halfHessian(hess, n, c, rx, ry, tx, ty, tdots, dists, L, la, ld); + + gsl_vector_free(x); + gsl_vector_free(y); + gsl_vector_free(rx); + gsl_vector_free(ry); + gsl_vector_free(tx); + gsl_vector_free(ty); + gsl_matrix_free(tdots); + gsl_matrix_free(dists); + gsl_permutation_free(indices_right); +} + + +void domain_energy_nakedHessian(gsl_matrix *hess, unsigned n, const gsl_vector *z, double c) { + + domain_energy_nakedHalfHessian(hess, n, z, c); + + #pragma omp parallel for + for (unsigned i = 1; i < 3 * n + 2; i++) { + for (unsigned j = 0; j < i; j++) { + gsl_matrix_set(hess, j, i, gsl_matrix_get(hess, i, j)); + } + } +} + + +double domain_energy_fixedEnergy(unsigned n, const gsl_vector *z, double c) { + double lagrangian; + + gsl_vector *x, *y; + + x = gsl_vector_alloc(n); + y = gsl_vector_alloc(n); + + #pragma omp parallel for + for (unsigned i = 0; i < n; i++) gsl_vector_set(x, i, gsl_vector_get(z, i)); + gsl_vector_set(y, 0, 0); + #pragma omp parallel for + for (unsigned i = 0; i < n - 1; i++) gsl_vector_set(y, i + 1, gsl_vector_get(z, i + n)); + double L = gsl_vector_get(z, 2 * n - 1); + + unsigned ii; + + gsl_vector *rx, *ry, *tx, *ty; + gsl_matrix *tdots, *dists; + + rx = gsl_vector_alloc(n); + ry = gsl_vector_alloc(n); + tx = gsl_vector_alloc(n); + ty = gsl_vector_alloc(n); + + tdots = gsl_matrix_alloc(n, n); + dists = gsl_matrix_alloc(n, n); + + domain_energy_rt(rx, n, x, 1); + domain_energy_rt(ry, n, y, 1); + domain_energy_rt(tx, n, x, -1); + domain_energy_rt(ty, n, y, -1); + + domain_energy_tdots(tdots, n, tx, ty); + domain_energy_dists(dists, n, rx, ry); + + lagrangian = domain_energy_energy(n, c, rx, ry, tx, ty, tdots, dists, L); + + gsl_vector_free(x); + gsl_vector_free(y); + gsl_vector_free(rx); + gsl_vector_free(ry); + gsl_vector_free(tx); + gsl_vector_free(ty); + + gsl_matrix_free(tdots); + gsl_matrix_free(dists); + + return lagrangian; +}; + + +// The fixed functions. + +double domain_energy_fixedLagrangian(unsigned n, const gsl_vector *z, double c) { + double lagrangian; + + gsl_vector *x, *y; + + x = gsl_vector_alloc(n); + y = gsl_vector_alloc(n); + + #pragma omp parallel for + for (unsigned i = 0; i < n; i++) gsl_vector_set(x, i, gsl_vector_get(z, i)); + gsl_vector_set(y, 0, 0); + #pragma omp parallel for + for (unsigned i = 0; i < n - 1; i++) gsl_vector_set(y, i + 1, gsl_vector_get(z, i + n)); + double L = gsl_vector_get(z, 2 * n - 1); + double la = gsl_vector_get(z, 2 * n); + double ld[n]; + #pragma omp parallel for + for (unsigned i = 0; i < n; i++) { + ld[i] = gsl_vector_get(z, 2 * n + 1 + i); + } + double lx = gsl_vector_get(z, 3 * n + 1); + double ly = gsl_vector_get(z, 3 * n + 2); + + unsigned ii; + + gsl_vector *rx, *ry, *tx, *ty; + gsl_matrix *tdots, *dists; + + rx = gsl_vector_alloc(n); + ry = gsl_vector_alloc(n); + tx = gsl_vector_alloc(n); + ty = gsl_vector_alloc(n); + + tdots = gsl_matrix_alloc(n, n); + dists = gsl_matrix_alloc(n, n); + + domain_energy_rt(rx, n, x, 1); + domain_energy_rt(ry, n, y, 1); + domain_energy_rt(tx, n, x, -1); + domain_energy_rt(ty, n, y, -1); + + domain_energy_tdots(tdots, n, tx, ty); + domain_energy_dists(dists, n, rx, ry); + + double area = domain_energy_area(n, x, y); + + lagrangian = domain_energy_lagrangian(n, c, rx, ry, tx, ty, tdots, dists, L, area, la, ld); + + double xtot = 0; + for (unsigned i = 0; i < n; i++) xtot += gsl_vector_get(z, i); + + double ytot = 0; + for (unsigned i = 1; i < n; i++) ytot += gsl_vector_get(z, i + n - 1); + + lagrangian += - lx * xtot - ly * ytot; + + gsl_vector_free(x); + gsl_vector_free(y); + gsl_vector_free(rx); + gsl_vector_free(ry); + gsl_vector_free(tx); + gsl_vector_free(ty); + + gsl_matrix_free(tdots); + gsl_matrix_free(dists); + + return lagrangian; +}; + + +void domain_energy_fixedGradient(gsl_vector *grad, unsigned n, const gsl_vector *z, + double c) { + + // Setting pointers to give the elements of z more convenient names. + double L = gsl_vector_get(z, 2 * n - 1); + double la = gsl_vector_get(z, 2 * n); + double ld[n]; + #pragma omp parallel for + for (unsigned i = 0; i < n; i++) { + ld[i] = gsl_vector_get(z, 2 * n + 1 + i); + } + double lx = gsl_vector_get(z, 3 * n + 1); + double ly = gsl_vector_get(z, 3 * n + 2); + + unsigned ii; + + gsl_vector *rx, *ry, *tx, *ty, *freegrad; + gsl_matrix *tdots, *dists; + + rx = gsl_vector_alloc(n); + ry = gsl_vector_alloc(n); + tx = gsl_vector_alloc(n); + ty = gsl_vector_alloc(n); + + tdots = gsl_matrix_alloc(n, n); + dists = gsl_matrix_alloc(n, n); + + freegrad = gsl_vector_alloc(3 * n + 2); + + double area = domain_energy_init(n, z, rx, ry, tx, ty, tdots, dists); + + domain_energy_gradient(freegrad, n, c, rx, ry, tx, ty, tdots, dists, L, area, la, ld); + + #pragma omp parallel for private(ii) + for (unsigned i =0; i < 3 * n + 2; i++) { + if (i != n) { + if (i < n) ii = i; + if (i > n) ii = i - 1; + + gsl_vector_set(grad, ii, gsl_vector_get(freegrad, i)); + } + } + + + #pragma omp parallel for private(ii) + for (unsigned i = 0; i < n; i++) { + gsl_vector_set(grad, i, gsl_vector_get(grad, i) - lx); + if (i != 0) gsl_vector_set(grad, n + i - 1, gsl_vector_get(grad, i + n - 1) -ly); + } + + + double xtot = 0; + for (unsigned i = 0; i < n; i++) xtot += gsl_vector_get(z, i); + double ytot = 0; + for (unsigned i = 1; i < n; i++) ytot += gsl_vector_get(z, i + n - 1); + + gsl_vector_set(grad, 3 * n + 1, -xtot); + + gsl_vector_set(grad, 3 * n + 2, -ytot); + + gsl_vector_free(rx); + gsl_vector_free(ry); + gsl_vector_free(tx); + gsl_vector_free(ty); + gsl_vector_free(freegrad); + gsl_matrix_free(tdots); + gsl_matrix_free(dists); + +} + + +void domain_energy_fixedHalfHessian(gsl_matrix *hess, unsigned n, const gsl_vector *z, + double c) { + + // Setting pointers to give the elements of z more convenient names. + double L = gsl_vector_get(z, 2 * n - 1); + double la = gsl_vector_get(z, 2 * n); + double ld[n]; + #pragma omp parallel for + for (unsigned i = 0; i < n; i++) { + ld[i] = gsl_vector_get(z, 2 * n + 1 + i); + } + double lx = gsl_vector_get(z, 3 * n + 1); + double ly = gsl_vector_get(z, 3 * n + 2); + + gsl_vector *rx, *ry, *tx, *ty; + gsl_matrix *tdots, *dists, *freehess; + + rx = gsl_vector_alloc(n); + ry = gsl_vector_alloc(n); + tx = gsl_vector_alloc(n); + ty = gsl_vector_alloc(n); + + tdots = gsl_matrix_alloc(n, n); + dists = gsl_matrix_alloc(n, n); + + freehess = gsl_matrix_alloc(3 * n + 2, 3 * n + 2); + + double area = domain_energy_init(n, z, rx, ry, tx, ty, tdots, dists); + + domain_energy_halfHessian(freehess, n, c, rx, ry, tx, ty, tdots, dists, L, la, ld); + + gsl_matrix *m_dxidxj, *m_dyidyj, *m_dxidxi, *m_dyidyi, *m_dxidxii, + *m_dyidyii, *m_dxidyj, *m_dxidyi, *m_dxiidyi, *m_dxidyii; + gsl_vector *v_ones, *v_storage; + gsl_permutation *indicesRight, *indicesLeft, *indices2Right; + + unsigned ii, jj, iii; + double rx_i, rx_j, tx_i, tx_j, tdt_ij, d_ij, rx_ii, rx_jj, tx_ii, tx_jj, + ry_i, ry_j, ty_i, ty_j, ry_ii, ry_jj, ty_ii, ty_jj, tdt_iij, d_iij, + tdt_ijj, d_ijj, tdt_iijj, d_iijj, rx_iii, tx_iii, ry_iii, ty_iii, + d_ij3, d_iij3, d_ijj3, d_iijj3, d_ij5, d_iij5, d_ijj5, d_iijj5; + + gsl_matrix_set_zero(hess); + + #pragma omp parallel for private(ii, jj) + for (unsigned i = 0; i < 3 * n + 2; i++) { + if (i != n) { + if (i < n) ii = i; + if (i > n) ii = i - 1; + for (unsigned j = 0; j <= i; j++) { + if (j != n) { + if (j < n) jj = j; + if (j > n) jj = j - 1; + + gsl_matrix_set(hess, ii, jj, gsl_matrix_get(freehess, i, j)); + } + } + } + } + + indicesRight = gsl_permutation_alloc(n); + indicesLeft = gsl_permutation_alloc(n); + indices2Right = gsl_permutation_alloc(n); + + gsl_permutation_init(indicesRight); + gsl_permutation_init(indicesLeft); + gsl_permutation_over(n, indicesRight, true); + gsl_permutation_over(n, indicesLeft, false); + gsl_permutation_memcpy(indices2Right, indicesRight); + gsl_permutation_over(n, indices2Right, true); + + + #pragma omp parallel for + for (unsigned i = 0; i < n; i++) { + gsl_matrix_set(hess, 3 * n + 1, i, -1); + if (i != 0) gsl_matrix_set(hess, 3 * n + 2, n + i - 1, -1); + } + + gsl_vector_free(rx); + gsl_vector_free(ry); + gsl_vector_free(tx); + gsl_vector_free(ty); + gsl_matrix_free(tdots); + gsl_matrix_free(dists); + gsl_matrix_free(freehess); + + gsl_permutation_free(indicesLeft); + gsl_permutation_free(indicesRight); + gsl_permutation_free(indices2Right); +} + + +void domain_energy_fixedHessian(gsl_matrix *hess, unsigned n, const gsl_vector *z, + double c) { + + domain_energy_fixedHalfHessian(hess, n, z, c); + + #pragma omp parallel for + for (unsigned i = 1; i < 3 * n + 3; i++) { + for (unsigned j = 0; j < i; j++) { + gsl_matrix_set(hess, j, i, gsl_matrix_get(hess, i, j)); + } + } +} + + +// The random functions. + +double domain_energy_randEnergy(unsigned n, const gsl_vector *z, + unsigned ord, const gsl_vector *k, const gsl_vector *a, + const gsl_vector *phi) { + + double energy; + + gsl_vector *x, *y; + + x = gsl_vector_alloc(n); + y = gsl_vector_alloc(n); + + #pragma omp parallel for + for (unsigned i = 0; i < n; i++) gsl_vector_set(x, i, gsl_vector_get(z, i)); + #pragma omp parallel for + for (unsigned i = 0; i < n; i++) gsl_vector_set(y, i, gsl_vector_get(z, i + n)); + + unsigned ii; + + gsl_vector *tx, *ty; + + tx = gsl_vector_alloc(n); + ty = gsl_vector_alloc(n); + + domain_energy_rt(tx, n, x, -1); + domain_energy_rt(ty, n, y, -1); + + energy = 0; + + for (unsigned i = 0; i < n; i++) { + for (unsigned j = 0; j < ord; j++) { + energy += gsl_vector_get(a, j) * gsl_sf_sin(gsl_vector_get(k, j) * gsl_vector_get(x, i) + gsl_vector_get(k, j + ord) * gsl_vector_get(y, i) + gsl_vector_get(phi, j)) * (gsl_vector_get(ty, i) / gsl_vector_get(k, j) - gsl_vector_get(tx, i) / gsl_vector_get(k, j + ord)) / 2; + } + } + gsl_vector_free(x); + gsl_vector_free(y); + gsl_vector_free(tx); + gsl_vector_free(ty); + + return energy; +}; + + +void domain_energy_randGradient(gsl_vector *grad, unsigned n, + const gsl_vector *z, unsigned ord, const gsl_vector *k, + const gsl_vector *a, const gsl_vector *phi) { + + gsl_permutation *indices_right; + indices_right = gsl_permutation_alloc(n); + + gsl_permutation_init(indices_right); + gsl_permutation_over(n, indices_right, true); + + gsl_vector *x, *y, *tx, *ty; + + x = gsl_vector_alloc(n); + y = gsl_vector_alloc(n); + + #pragma omp parallel for + for (unsigned i = 0; i < n; i++) gsl_vector_set(x, i, gsl_vector_get(z, i)); + #pragma omp parallel for + for (unsigned i = 0; i < n; i++) gsl_vector_set(y, i, gsl_vector_get(z, i + n)); + + unsigned ii; + + tx = gsl_vector_alloc(n); + ty = gsl_vector_alloc(n); + + domain_energy_rt(tx, n, x, -1); + domain_energy_rt(ty, n, y, -1); + + double thesumx, thesumy, aj, kxj, kyj, phij, xi, yi, xii, yii, txi, tyi; + + #pragma omp parallel for private(ii, xi, yi, xii, yii, txi, tyi, thesumx, thesumy, aj, kxj, kyj, phij) + for (unsigned i = 0; i < n; i++) { + ii = gsl_permutation_get(indices_right, i); + + xi = gsl_vector_get(x, i); + yi = gsl_vector_get(y, i); + xii = gsl_vector_get(x, ii); + yii = gsl_vector_get(y, ii); + txi = gsl_vector_get(tx, i); + tyi = gsl_vector_get(ty, i); + + thesumx = 0; + thesumy = 0; + + for (unsigned j = 0; j < ord; j++) { + aj = gsl_vector_get(a, j); + kxj = gsl_vector_get(k, j); + kyj = gsl_vector_get(k, ord + j); + phij = gsl_vector_get(phi, j); + + thesumx += aj * (kxj * gsl_sf_cos(kxj * xi + kyj * yi + phij) * (tyi / kxj - txi / kyj) + + (gsl_sf_sin(kxj * xi + kyj * yi + phij) - gsl_sf_sin(kxj * xii + kyj * yii + phij)) / kyj); + + thesumy += aj * (kyj * gsl_sf_cos(kxj * xi + kyj * yi + phij) * (tyi / kxj - txi / kyj) - + (gsl_sf_sin(kxj * xi + kyj * yi + phij) - gsl_sf_sin(kxj * xii + kyj * yii + phij)) / kxj); + } + + gsl_vector_set(grad, i, gsl_vector_get(grad, i) + thesumx / 2); + gsl_vector_set(grad, n + i, gsl_vector_get(grad, n + i) + thesumy / 2); + } + + gsl_vector_free(x); + gsl_vector_free(y); + gsl_vector_free(tx); + gsl_vector_free(ty); + gsl_permutation_free(indices_right); +} + + +void domain_energy_randHalfHessian(gsl_matrix *hess, unsigned n, + const gsl_vector *z, unsigned ord, const gsl_vector *k, const gsl_vector *a, + const gsl_vector *phi) { + + gsl_permutation *indices_right; + indices_right = gsl_permutation_alloc(n); + + gsl_permutation_init(indices_right); + gsl_permutation_over(n, indices_right, true); + + gsl_vector *x, *y, *tx, *ty; + + x = gsl_vector_alloc(n); + y = gsl_vector_alloc(n); + + #pragma omp parallel for + for (unsigned i = 0; i < n; i++) gsl_vector_set(x, i, gsl_vector_get(z, i)); + #pragma omp parallel for + for (unsigned i = 0; i < n; i++) gsl_vector_set(y, i, gsl_vector_get(z, i + n)); + + unsigned ii; + + tx = gsl_vector_alloc(n); + ty = gsl_vector_alloc(n); + + domain_energy_rt(tx, n, x, -1); + domain_energy_rt(ty, n, y, -1); + + double thesumx, thesumy, thesumxy, thesumxx, thesumyy, thesumxxy, thesumxyy, aj, kxj, kyj, phij, xi, yi, xii, yii, txi, tyi; + + #pragma omp parallel for private(ii, xi, yi, xii, yii, txi, tyi, thesumx, thesumy, thesumxy, thesumxx, thesumyy, thesumxxy, thesumxyy, aj, kxj, kyj, phij) + for (unsigned i = 0; i < n; i++) { + ii = gsl_permutation_get(indices_right, i); + + xi = gsl_vector_get(x, i); + yi = gsl_vector_get(y, i); + xii = gsl_vector_get(x, ii); + yii = gsl_vector_get(y, ii); + txi = gsl_vector_get(tx, i); + tyi = gsl_vector_get(ty, i); + + thesumx = 0; + thesumy = 0; + thesumxy = 0; + thesumxx = 0; + thesumyy = 0; + thesumxxy = 0; + thesumxyy = 0; + + for (unsigned j = 0; j < ord; j++) { + aj = gsl_vector_get(a, j); + kxj = gsl_vector_get(k, j); + kyj = gsl_vector_get(k, ord + j); + phij = gsl_vector_get(phi, j); + + thesumx += aj * ( - gsl_pow_2(kxj) * gsl_sf_sin(kxj * xi + kyj * yi + phij) * (tyi / kxj - txi / kyj) + + 2 * gsl_sf_cos(kxj * xi + kyj * yi + phij) * kxj / kyj); + + thesumy += aj * ( - gsl_pow_2(kyj) * gsl_sf_sin(kxj * xi + kyj * yi + phij) * (tyi / kxj - txi / kyj) - + 2 * gsl_sf_cos(kxj * xi + kyj * yi + phij) * kyj / kxj); + + thesumxy += - aj * kxj * kyj * gsl_sf_sin(kxj * xi + kyj * yi + phij) * (tyi / kxj - txi / kyj); + + thesumxx += - aj * gsl_sf_cos(kxj * xii + kyj * yii + phij) * kxj / kyj; + + thesumyy += aj * gsl_sf_cos(kxj * xii + kyj * yii + phij) * kyj / kxj; + + thesumxyy += - aj * gsl_sf_cos(kxj * xii + kyj * yii + phij); + + thesumxxy += aj * gsl_sf_cos(kxj * xii + kyj * yii + phij); + } + + gsl_matrix_set(hess, i, i, gsl_matrix_get(hess, i, i) + thesumx / 2); + gsl_matrix_set(hess, n + i, n + i, gsl_matrix_get(hess, n + i, n + i) + thesumy / 2); + gsl_matrix_set(hess, i + n, i, gsl_matrix_get(hess, n + i, i) + thesumxy / 2); + if (i == 0) { + gsl_matrix_set(hess, n - 1, 0, gsl_matrix_get(hess, n - 1, 0) + thesumxx / 2); + gsl_matrix_set(hess, 2 * n - 1, n, gsl_matrix_get(hess, 2 * n - 1, n) + thesumyy / 2); + } else { + gsl_matrix_set(hess, i, ii, gsl_matrix_get(hess, i, ii) + thesumxx / 2); + gsl_matrix_set(hess, n + i, n + ii, gsl_matrix_get(hess, n + i, n + ii) + thesumyy / 2); + } + gsl_matrix_set(hess, n + i, ii, gsl_matrix_get(hess, n + i, ii) + thesumxxy / 2); + gsl_matrix_set(hess, n + ii, i, gsl_matrix_get(hess, n + ii, i) + thesumxyy / 2); + } + gsl_vector_free(x); + gsl_vector_free(y); + gsl_vector_free(tx); + gsl_vector_free(ty); + gsl_permutation_free(indices_right); +} + + +// The random naked functions. + +double domain_energy_nakedRandLagrangian(unsigned n, const gsl_vector *z, + double c, unsigned ord, const gsl_vector *k, const gsl_vector *a, + const gsl_vector *phi) { + + double lagrangian, randEnergy; + + lagrangian = domain_energy_nakedLagrangian(n, z, c); + randEnergy = domain_energy_randEnergy(n, z, ord, k, a, phi); + + return lagrangian + randEnergy; +} + + +void domain_energy_nakedRandGradient(gsl_vector *grad, unsigned n, + const gsl_vector *z, double c, unsigned ord, const gsl_vector *k, + const gsl_vector *a, const gsl_vector *phi) { + + domain_energy_nakedGradient(grad, n, z, c); + domain_energy_randGradient(grad, n, z, ord, k, a, phi); +} + + +void domain_energy_nakedRandHessian(gsl_matrix *hess, unsigned n, + const gsl_vector *z, double c, unsigned ord, const gsl_vector *k, + const gsl_vector *a, const gsl_vector *phi) { + + domain_energy_nakedHalfHessian(hess, n, z, c); + domain_energy_randHalfHessian(hess, n, z, ord, k, a, phi); + + #pragma omp parallel for + for (unsigned i = 1; i < 3 * n + 2; i++) { + for (unsigned j = 0; j < i; j++) { + gsl_matrix_set(hess, j, i, gsl_matrix_get(hess, i, j)); + } + } +} + + +// The well functions. + +double domain_energy_wellEnergy(unsigned n, const gsl_vector *z, double w, + double s) { + + double xi, yi, txi, tyi, energy; + gsl_vector *x, *y, *tx, *ty; + + x = gsl_vector_alloc(n); + y = gsl_vector_alloc(n); + tx = gsl_vector_alloc(n); + ty = gsl_vector_alloc(n); + + #pragma omp parallel for + for (unsigned i = 0; i < n; i++) gsl_vector_set(x, i, gsl_vector_get(z, i)); + #pragma omp parallel for + for (unsigned i = 0; i < n; i++) gsl_vector_set(y, i, + gsl_vector_get(z, i + n)); + + domain_energy_rt(tx, n, x, -1); + domain_energy_rt(ty, n, y, -1); + + energy = 0; + + for (unsigned i = 0; i < n; i++) { + xi = gsl_vector_get(x, i); + yi = gsl_vector_get(y, i); + txi = gsl_vector_get(tx, i); + tyi = gsl_vector_get(ty, i); + + energy += ((gsl_sf_exp(s * (xi - w)) - gsl_sf_exp(-s * (xi + w))) * tyi - + (gsl_sf_exp(s * (yi - w)) - gsl_sf_exp(-s * (yi + w))) * txi) / s; + } + + gsl_vector_free(x); + gsl_vector_free(y); + gsl_vector_free(tx); + gsl_vector_free(ty); + + return energy; +}; + + +void domain_energy_wellGradient(gsl_vector *grad, unsigned n, + const gsl_vector *z, double w, double s) { + + unsigned ii; + double xi, yi, xii, yii, txi, tyi, grad_xi, grad_yi; + gsl_vector *x, *y, *tx, *ty; + gsl_permutation *indices_right; + + x = gsl_vector_alloc(n); + y = gsl_vector_alloc(n); + tx = gsl_vector_alloc(n); + ty = gsl_vector_alloc(n); + indices_right = gsl_permutation_alloc(n); + + #pragma omp parallel for + for (unsigned i = 0; i < n; i++) gsl_vector_set(x, i, gsl_vector_get(z, i)); + #pragma omp parallel for + for (unsigned i = 0; i < n; i++) gsl_vector_set(y, i, gsl_vector_get(z, i + n)); + + gsl_permutation_init(indices_right); + gsl_permutation_over(n, indices_right, true); + + domain_energy_rt(tx, n, x, -1); + domain_energy_rt(ty, n, y, -1); + + + #pragma omp parallel for private(ii, xi, yi, xii, yii, txi, tyi, grad_xi,\ + grad_yi) + for (unsigned i = 0; i < n; i++) { + ii = gsl_permutation_get(indices_right, i); + + xi = gsl_vector_get(x, i); + yi = gsl_vector_get(y, i); + xii = gsl_vector_get(x, ii); + yii = gsl_vector_get(y, ii); + txi = gsl_vector_get(tx, i); + tyi = gsl_vector_get(ty, i); + + grad_xi = tyi * (gsl_sf_exp(s * (xi - w)) + gsl_sf_exp(-s * (xi + w))) + + (gsl_sf_exp(s * (yi - w)) - gsl_sf_exp(-s * (yi + w)) - + gsl_sf_exp(s * (yii - w)) + gsl_sf_exp(-s * (yii + w))) / s; + + grad_yi = - txi * (gsl_sf_exp(s * (yi - w)) + gsl_sf_exp(-s * (yi + w))) + + (- gsl_sf_exp(s * (xi - w)) + gsl_sf_exp(-s * (xi + w)) + + gsl_sf_exp(s * (xii - w)) - gsl_sf_exp(-s * (xii + w))) / s; + + gsl_vector_set(grad, i, gsl_vector_get(grad, i) + grad_xi); + + gsl_vector_set(grad, n + i, gsl_vector_get(grad, n + i) + grad_yi); + } + + gsl_vector_free(x); + gsl_vector_free(y); + gsl_vector_free(tx); + gsl_vector_free(ty); + gsl_permutation_free(indices_right); +} + + +void domain_energy_wellHalfHessian(gsl_matrix *hess, unsigned n, + const gsl_vector *z, double w, double s) { + + unsigned ii; + double xi, yi, xii, yii, txi, tyi, hess_xi, hess_yi, hess_xiyi, hess_xiiyi, + hess_xiyii, exp_mxi, exp_pxi, exp_myi, exp_pyi, exp_mxii, exp_myii, + exp_pxii, exp_pyii; + gsl_vector *x, *y, *tx, *ty; + gsl_permutation *indices_right; + + x = gsl_vector_alloc(n); + y = gsl_vector_alloc(n); + tx = gsl_vector_alloc(n); + ty = gsl_vector_alloc(n); + indices_right = gsl_permutation_alloc(n); + + for (unsigned i = 0; i < n; i++) gsl_vector_set(x, i, gsl_vector_get(z, i)); + for (unsigned i = 0; i < n; i++) { + gsl_vector_set(y, i, gsl_vector_get(z, i + n)); + } + + domain_energy_rt(tx, n, x, -1); + domain_energy_rt(ty, n, y, -1); + + gsl_permutation_init(indices_right); + gsl_permutation_over(n, indices_right, true); + + #pragma omp parallel for private(ii, xi, yi, xii, yii, txi, tyi, hess_xi,\ + hess_yi, hess_xiyi, hess_xiiyi, hess_xiyii, exp_mxi, exp_pxi, exp_myi,\ + exp_pyi, exp_mxii, exp_myii, exp_pxii, exp_pyii) + for (unsigned i = 0; i < n; i++) { + ii = gsl_permutation_get(indices_right, i); + + xi = gsl_vector_get(x, i); + yi = gsl_vector_get(y, i); + xii = gsl_vector_get(x, ii); + yii = gsl_vector_get(y, ii); + txi = gsl_vector_get(tx, i); + tyi = gsl_vector_get(ty, i); + exp_mxi = gsl_sf_exp(-s * (xi + w)); + exp_pxi = gsl_sf_exp(s * (xi - w)); + exp_myi = gsl_sf_exp(-s * (yi + w)); + exp_pyi = gsl_sf_exp(s * (yi - w)); + exp_mxii = gsl_sf_exp(-s * (xii + w)); + exp_pxii = gsl_sf_exp(s * (xii - w)); + exp_myii = gsl_sf_exp(-s * (yii + w)); + exp_pyii = gsl_sf_exp(s * (yii - w)); + + hess_xi = tyi * s * (exp_pxi - exp_mxi); + hess_yi = - txi * s * (exp_pyi - exp_myi); + + hess_xiyi = - exp_pxi - exp_mxi + exp_pyi + exp_myi; + + hess_xiyii = - exp_pyii - exp_myii; + hess_xiiyi = exp_pxii + exp_mxii; + + gsl_matrix_set(hess, i, i, gsl_matrix_get(hess, i, i) + hess_xi); + gsl_matrix_set(hess, n + i, n + i, gsl_matrix_get(hess, n + i, n + i) + + hess_yi); + gsl_matrix_set(hess, i + n, i, gsl_matrix_get(hess, n + i, i) + hess_xiyi); + gsl_matrix_set(hess, n + i, ii, gsl_matrix_get(hess, n + i, ii) + + hess_xiiyi); + gsl_matrix_set(hess, n + ii, i, gsl_matrix_get(hess, n + ii, i) + + hess_xiyii); + } + + gsl_vector_free(x); + gsl_vector_free(y); + gsl_vector_free(tx); + gsl_vector_free(ty); + gsl_permutation_free(indices_right); +} + + +// The naked well functions. + +double domain_energy_nakedWellLagrangian(unsigned n, const gsl_vector *z, + double c, double w, double s) { + + double nakedLagrangian, wellEnergy; + + nakedLagrangian = domain_energy_nakedLagrangian(n, z, c); + wellEnergy = domain_energy_wellEnergy(n, z, w, s); + + return nakedLagrangian + wellEnergy; +} + + +void domain_energy_nakedWellGradient(gsl_vector *grad, unsigned n, + const gsl_vector *z, double c, double w, double s) { + + domain_energy_nakedGradient(grad, n, z, c); + domain_energy_wellGradient(grad, n, z, w, s); +} + + +void domain_energy_nakedWellHessian(gsl_matrix *hess, unsigned n, + const gsl_vector *z, double c, double w, double s) { + + domain_energy_nakedHalfHessian(hess, n, z, c); + domain_energy_wellHalfHessian(hess, n, z, w, s); + + #pragma omp parallel for + for (unsigned i = 1; i < 3 * n + 2; i++) { + for (unsigned j = 0; j < i; j++) { + gsl_matrix_set(hess, j, i, gsl_matrix_get(hess, i, j)); + } + } +} + + +// The random well functions. + +double domain_energy_randWellLagrangian(unsigned n, const gsl_vector *z, + double c, unsigned ord, const gsl_vector *k, const gsl_vector *a, + const gsl_vector *phi, double w, double s) { + + double lagrangian, randEnergy, wellEnergy; + + lagrangian = domain_energy_nakedLagrangian(n, z, c); + randEnergy = domain_energy_randEnergy(n, z, ord, k, a, phi); + wellEnergy = domain_energy_wellEnergy(n, z, w, s); + + return lagrangian + randEnergy + wellEnergy; +} + + +void domain_energy_randWellGradient(gsl_vector *grad, unsigned n, + const gsl_vector *z, double c, unsigned ord, const gsl_vector *k, + const gsl_vector *a, const gsl_vector *phi, double w, double s) { + + domain_energy_nakedGradient(grad, n, z, c); + domain_energy_randGradient(grad, n, z, ord, k, a, phi); + domain_energy_wellGradient(grad, n, z, w, s); +} + + +void domain_energy_randWellHessian(gsl_matrix *hess, unsigned n, + const gsl_vector *z, double c, unsigned ord, const gsl_vector *k, + const gsl_vector *a, const gsl_vector *phi, double w, double s) { + + domain_energy_nakedHalfHessian(hess, n, z, c); + domain_energy_randHalfHessian(hess, n, z, ord, k, a, phi); + domain_energy_wellHalfHessian(hess, n, z, w, s); + + #pragma omp parallel for + for (unsigned i = 1; i < 3 * n + 2; i++) { + for (unsigned j = 0; j < i; j++) { + gsl_matrix_set(hess, j, i, gsl_matrix_get(hess, i, j)); + } + } +} diff --git a/src/domain_energy.h b/src/domain_energy.h new file mode 100644 index 0000000..6cdf882 --- /dev/null +++ b/src/domain_energy.h @@ -0,0 +1,41 @@ +#ifndef DOMAIN_ENERGY_H +#define DOMAIN_ENERGY_H + +#include <gsl/gsl_vector.h> +#include <gsl/gsl_matrix.h> + +double domain_energy_nakedEnergy(unsigned n, const gsl_vector *z, double c); + +double domain_energy_nakedLagrangian(unsigned n, const gsl_vector *z, double c); + +double domain_energy_nakedGradient(gsl_vector *grad, unsigned n, const gsl_vector *z, double c); + +double domain_energy_nakedHessian(gsl_matrix *hess, unsigned n, const gsl_vector *z, double c); + +double domain_energy_fixedEnergy(unsigned n, const gsl_vector *z, double c); + +double domain_energy_fixedLagrangian(unsigned n, const gsl_vector *z, double c); + +double domain_energy_fixedGradient(gsl_vector *grad, unsigned n, const gsl_vector *z, double c); + +double domain_energy_fixedHessian(gsl_matrix *hess, unsigned n, const gsl_vector *z, double c); + +double domain_energy_nakedRandLagrangian(unsigned n, const gsl_vector *z, double c, unsigned ord, const gsl_vector *k, const gsl_vector *a, const gsl_vector *phi); + +double domain_energy_nakedRandGradient(gsl_vector *grad, unsigned n, const gsl_vector *z, double c, unsigned ord, const gsl_vector *k, const gsl_vector *a, const gsl_vector *phi); + +double domain_energy_nakedRandHessian(gsl_matrix *hess, unsigned n, const gsl_vector *z, double c, unsigned ord, const gsl_vector *k, const gsl_vector *a, const gsl_vector *phi); + +double domain_energy_nakedWellLagrangian(unsigned n, const gsl_vector *z, double c, double w, double s); + +double domain_energy_nakedWellGradient(gsl_vector *grad, unsigned n, const gsl_vector *z, double c, double w, double s); + +double domain_energy_nakedWellHessian(gsl_matrix *hess, unsigned n, const gsl_vector *z, double c, double w, double s); + +double domain_energy_randWellLagrangian(unsigned n, const gsl_vector *z, double c, unsigned ord, const gsl_vector *k, const gsl_vector *a, const gsl_vector *phi, double w, double s); + +double domain_energy_randWellGradient(gsl_vector *grad, unsigned n, const gsl_vector *z, double c, unsigned ord, const gsl_vector *k, const gsl_vector *a, const gsl_vector *phi, double w, double s); + +double domain_energy_randWellHessian(gsl_matrix *hess, unsigned n, const gsl_vector *z, double c, unsigned ord, const gsl_vector *k, const gsl_vector *a, const gsl_vector *phi, double w, double s); + +#endif diff --git a/src/domain_improve.cpp b/src/domain_improve.cpp new file mode 100644 index 0000000..f69e85f --- /dev/null +++ b/src/domain_improve.cpp @@ -0,0 +1,145 @@ +/* domain_improve.cpp + * + * Copyright (C) 2013 Jaron Kent-Dobias + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +/* A program which facilitates automated mapping of bifurcation points in the + * energy of a system where the Hessian is available. Currently, only a one + * dimensional parameter space is supported. + */ + +#include "domain_energy.h" +#include "domain_minimize.h" +#include "domain_eigen.h" + +#include <unistd.h> +#include <stdio.h> +#include <iostream> +#include <stdlib.h> +#include <math.h> +#include <string> + +// GSL includes. +#include <gsl/gsl_sf.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_complex.h> +#include <gsl/gsl_complex_math.h> +#include <gsl/gsl_vector.h> +#include <gsl/gsl_permutation.h> +#include <gsl/gsl_permute_vector.h> +#include <gsl/gsl_blas.h> +#include <gsl/gsl_linalg.h> +#include <gsl/gsl_sort_vector.h> + + +// Initializes the program. +int main(int argc, char *argv[]) { + + int opt, min_fails; + unsigned n, N, num; + double c, g0, g, eps, energy; + char *filename; + bool eigenpres = true; + + // Setting default values. + eps = 0; + num = 25; + + gsl_vector *z, *old_z, *eigenvalues, *trueEigenvalues; + gsl_permutation *eigenorder, *trueEigenorder; + + while ((opt = getopt(argc, argv, "n:c:d:g:h:i:N:p:m:j:e:t:s")) != -1) { + switch (opt) { + case 'n': + n = atoi(optarg); + break; + case 'N': + N = atoi(optarg); + break; + case 'g': + g0 = atof(optarg); + break; + case 'i': + filename = optarg; + break; + case 'e': + eps = atof(optarg); + break; + default: + exit(EXIT_FAILURE); + } + } + + z = gsl_vector_alloc(3 * n + 3); + old_z = gsl_vector_alloc(3 * n + 3); + eigenvalues = gsl_vector_alloc(3 * n + 3); + trueEigenvalues = gsl_vector_alloc(3 * n + 4); + eigenorder = gsl_permutation_alloc(3 * n + 3); + trueEigenorder = gsl_permutation_alloc(3 * n + 4); + + g = g0; + + char ch; + double throwaway; + + FILE *f = fopen(filename, "r+"); + while (ch != '\n') ch = fgetc(f); + ch = 'a'; + while (ch != '\n' && ch != '\t') ch = fgetc(f); + if (ch == '\n') eigenpres = false; + + rewind(f); + + fscanf(f, "%le\t", &c); + fscanf(f, "%le\n", &energy); + + if (eigenpres) { + ch = 'a'; + while (ch != '\n') ch = fgetc(f); + } + gsl_vector_fscanf(f, z); + fclose(f); + + min_fails = domain_minimize_fixed(z, n, c, eps, g, N, 4, 2); + + if (min_fails) { + printf("BIFUR: Initial relaxation failed, exiting.\n"); + return 1; + } + +// domain_eigen_values(eigenvalues, n, z, c); +// domain_eigen_sort(eigenorder, n, num, eigenvalues); + + energy = domain_energy_fixedEnergy(n, z, c); + unsigned ii; + + FILE *newf = fopen(filename, "w"); + fprintf(newf, "%.12le\t%.12le\n", c, energy); + for (unsigned i = 0; i < num; i++) { + ii = gsl_permutation_get(eigenorder, i); + fprintf(newf, "%.12le\t", gsl_vector_get(eigenvalues, ii)); + } + fprintf(newf, "\n"); + for (unsigned i = 0; i < num; i++) { + ii = gsl_permutation_get(trueEigenorder, i); + fprintf(newf, "%.12le\t", gsl_vector_get(trueEigenvalues, ii)); + } + fprintf(newf, "\n"); + for (unsigned i = 0; i < 3 * n + 3; i++) { + fprintf(newf, "%.12le\t", gsl_vector_get(z, i)); + } + fclose(newf); +} diff --git a/src/domain_increase.cpp b/src/domain_increase.cpp new file mode 100644 index 0000000..9c95d8d --- /dev/null +++ b/src/domain_increase.cpp @@ -0,0 +1,120 @@ +/* domain_init.cpp + * + * Copyright (C) 2013 Jaron Kent-Dobias + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +// Initializes a circular domain, or a circular domain with a perturbation. + +#include <unistd.h> +#include <iostream> +#include <string> +#include <stdlib.h> +#include <math.h> + +// GSL includes. +#include <gsl/gsl_sf.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_vector.h> + +#include "domain_minimize.h" +#include "domain_energy.h" + + + + +int main(int argc, char *argv[]) { + + // Declaring variables. + gsl_vector *z, *new_z; + int opt; + unsigned n; + double c; + char out_file[20]; + + // Default values. + sprintf(out_file, "%s", "out.dat"); + + // GNU getopt in action. + while ((opt = getopt(argc, argv, "n:c:i:")) != -1) { + switch (opt) { + case 'n': + n = atoi(optarg); + break; + case 'c': + c = atof(optarg); + break; + case 'i': + sprintf(out_file, "%s", optarg); + break; + default: + exit(EXIT_FAILURE); + } + } + + z = gsl_vector_alloc(3 * n + 3); + new_z = gsl_vector_alloc(6 * n + 3); + + FILE *f = fopen(out_file, "r"); + gsl_vector_fscanf(f, z); + fclose(f); + + gsl_vector_set(new_z, 2 * n - 1, (gsl_vector_get(z, 0) + gsl_vector_get(z, n - 1)) / 2); + gsl_vector_set(new_z, 2 * n - 2, gsl_vector_get(z, n - 2)); + + for (unsigned i = 0; i < n - 1; i ++) { + gsl_vector_set(new_z, 2 * i, gsl_vector_get(z, i)); + gsl_vector_set(new_z, 2 * i + 1, + (gsl_vector_get(z, i) + gsl_vector_get(z, i + 1)) / 2); + } + + gsl_vector_set(new_z, 2 * n, gsl_vector_get(z, n) / 2); + gsl_vector_set(new_z, 2 * n + 1, gsl_vector_get(z, n)); + gsl_vector_set(new_z, 4 * n - 2, gsl_vector_get(z, 2 * n - 2) / 2); + gsl_vector_set(new_z, 4 * n - 3, gsl_vector_get(z, 2 * n - 2)); + + for (unsigned i = 1; i < n - 1; i ++) { + gsl_vector_set(new_z, 2 * (n + i) + 1, gsl_vector_get(z, n + i)); + gsl_vector_set(new_z, 2 * (n + i) , + (gsl_vector_get(z, n + i - 1) + gsl_vector_get(z, n + i)) / 2); + } + + gsl_vector_set(new_z, 4 * n - 1, gsl_vector_get(z, 2 * n -1)); + + gsl_vector_set(new_z, 4 * n, gsl_vector_get(z, 2 * n)); + + for (unsigned i = 0; i < n; i++) { + gsl_vector_set(new_z, 4 * n + 1 + 2 * i, gsl_vector_get(z, 2 * n + 1 + i)); + gsl_vector_set(new_z, 4 * n + 2 + 2 * i, gsl_vector_get(z, 2 * n + 1 + i)); + } + + gsl_vector_set(new_z, 6 * n + 1, gsl_vector_get(z, 3 * n + 1)); + gsl_vector_set(new_z, 6 * n + 2, gsl_vector_get(z, 3 * n + 2)); + + int result; + result = domain_minimize_fixed(new_z, 2 * n, c, 1e-8, 0.00001, 5000, 4, 2); + + if (result) printf("Converging failed."); + else { + FILE *fout = fopen(out_file, "w"); + gsl_vector_fprintf(fout, new_z, "%.10g"); + fclose(fout); + } + + gsl_vector_free(z); + gsl_vector_free(new_z); + +} + diff --git a/src/domain_minimize.cpp b/src/domain_minimize.cpp new file mode 100644 index 0000000..f316e4b --- /dev/null +++ b/src/domain_minimize.cpp @@ -0,0 +1,300 @@ +/* domain_minimize.cpp + * + * Copyright (C) 2013 Jaron Kent-Dobias + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +// A Newton's method solver for modulated domains. + +// GSL includes. +#include <gsl/gsl_math.h> +#include <gsl/gsl_vector.h> +#include <gsl/gsl_blas.h> + +// Gives the necessary functions for the Lagrangian, gradient, and Hessian. +#include "domain_energy.h" +#include "domain_newton.h" + +struct nakedgetgrad { + nakedgetgrad(double c, unsigned n): c(c), n(n) {} + void operator()(gsl_vector* grad, unsigned size, gsl_vector *state) {domain_energy_nakedGradient(grad, n, state, c);} + + private: + double c; + unsigned n; +}; + +struct nakedgethess { + nakedgethess(double c, unsigned n): c(c), n(n) {} + void operator()(gsl_matrix* hess, unsigned size, gsl_vector *state) {domain_energy_nakedHessian(hess, n, state, c);} + + private: + double c; + unsigned n; +}; + +struct nakedgetenergy { + nakedgetenergy(double c, unsigned n): c(c), n(n) {} + double operator()(unsigned size, gsl_vector *state) {return domain_energy_nakedLagrangian(n, state, c);} + + private: + double c; + unsigned n; +}; + +// Carries out Newton's method. +int domain_minimize_naked(gsl_vector *z, unsigned n, double c, double eps, unsigned N, double beta, double s, double sigma, double gamma, double eta0, bool verb) { + + unsigned size = 3 * n + 2; + unsigned params = 2 * n + 1; + nakedgetgrad grad(c, n); + nakedgethess hess(c, n); + nakedgetenergy energy(c, n); + + return domain_newton(z, size, params, energy, grad, hess, eps, N, beta, s, sigma, gamma, eta0, 0.1, 100, verb); +} + +struct fixedgetgrad { + fixedgetgrad(double c, unsigned n): c(c), n(n) {} + void operator()(gsl_vector* grad, unsigned size, gsl_vector *state) {domain_energy_fixedGradient(grad, n, state, c);} + + private: + double c; + unsigned n; +}; + +struct fixedgethess { + fixedgethess(double c, unsigned n): c(c), n(n) {} + void operator()(gsl_matrix* hess, unsigned size, gsl_vector *state) {domain_energy_fixedHessian(hess, n, state, c);} + + private: + double c; + unsigned n; +}; + +struct fixedgetenergy { + fixedgetenergy(double c, unsigned n): c(c), n(n) {} + double operator()(unsigned size, gsl_vector *state) {return domain_energy_fixedLagrangian(n, state, c);} + + private: + double c; + unsigned n; +}; + +// Carries out Newton's method. +int domain_minimize_fixed(gsl_vector *z, unsigned n, double c, double eps, unsigned N, double beta, double s, double sigma) { + + unsigned size = 3 * n + 3; + unsigned params = 2 * n; + fixedgetgrad grad(c, n); + fixedgethess hess(c, n); + fixedgetenergy energy(c, n); + + return domain_newton(z, size, params, energy, grad, hess, eps, N, beta, s, sigma, 0, 0, 0.1, 10, true); +} + +struct randgetgrad { + randgetgrad(double c, unsigned n, unsigned ord, const gsl_vector *k, const gsl_vector *a, const gsl_vector *phi): c(c), n(n), ord(ord), k(k), a(a), phi(phi) {} + void operator()(gsl_vector* grad, unsigned size, gsl_vector *state) {domain_energy_nakedRandGradient(grad, n, state, c, ord, k, a, phi);} + + private: + double c; + unsigned n; + unsigned ord; + const gsl_vector *k; + const gsl_vector *a; + const gsl_vector *phi; +}; + +struct randgethess { + randgethess(double c, unsigned n, unsigned ord, const gsl_vector *k, const gsl_vector *a, const gsl_vector *phi): c(c), n(n), ord(ord), k(k), a(a), phi(phi) {} + void operator()(gsl_matrix* hess, unsigned size, gsl_vector *state) {domain_energy_nakedRandHessian(hess, n, state, c, ord, k, a, phi);} + + private: + double c; + unsigned n; + unsigned ord; + const gsl_vector *k; + const gsl_vector *a; + const gsl_vector *phi; +}; + +struct randgetenergy { + randgetenergy(double c, unsigned n, unsigned ord, const gsl_vector *k, const gsl_vector *a, const gsl_vector *phi): c(c), n(n), ord(ord), k(k), a(a), phi(phi) {} + double operator()(unsigned size, gsl_vector *state) {return domain_energy_nakedRandLagrangian(n, state, c, ord, k, a, phi);} + + private: + double c; + unsigned n; + unsigned ord; + const gsl_vector *k; + const gsl_vector *a; + const gsl_vector *phi; +}; + +// Carries out Newton's method. +int domain_minimize_rand(gsl_vector *z, unsigned n, double c, unsigned ord, const gsl_vector *k, const gsl_vector *a, const gsl_vector *phi, double eps, unsigned N, double beta, double s, double sigma, double gamma, double bound, bool verb) { + + unsigned size = 3 * n + 2; + unsigned params = 2 * n + 1; + randgetgrad grad(c, n, ord, k, a, phi); + randgethess hess(c, n, ord, k, a, phi); + randgetenergy energy(c, n, ord, k, a, phi); + + return domain_newton(z, size, params, energy, grad, hess, eps, N, beta, s, sigma, gamma, bound, 0.1, 2, verb); +} + +struct nakedwellgetgrad { + nakedwellgetgrad(double c, unsigned n, double w, double s): c(c), n(n), w(w), s(s) {} + void operator()(gsl_vector* grad, unsigned size, gsl_vector *state) {domain_energy_nakedWellGradient(grad, n, state, c, w, s);} + + private: + double c; + double s; + double w; + unsigned n; +}; + +struct nakedwellgethess { + nakedwellgethess(double c, unsigned n, double w, double s): c(c), n(n), w(w), s(s) {} + void operator()(gsl_matrix* hess, unsigned size, gsl_vector *state) {domain_energy_nakedWellHessian(hess, n, state, c, w, s);} + + private: + double c; + double s; + double w; + unsigned n; +}; + +struct nakedwellgetenergy { + nakedwellgetenergy(double c, unsigned n, double w, double s): c(c), n(n), w(w), s(s) {} + double operator()(unsigned size, gsl_vector *state) {return domain_energy_nakedWellLagrangian(n, state, c, w, s);} + + private: + double c; + double s; + double w; + unsigned n; +}; + +// Carries out Newton's method. +int domain_minimize_nakedWell(gsl_vector *z, unsigned n, double c, double w, double ss, double eps, unsigned N, double beta, double s, double sigma, double gamma, double eta0, bool verb) { + + unsigned size = 3 * n + 2; + unsigned params = 2 * n + 1; + nakedwellgetgrad grad(c, n, w, ss); + nakedwellgethess hess(c, n, w, ss); + nakedwellgetenergy energy(c, n, w, ss); + + return domain_newton(z, size, params, energy, grad, hess, eps, N, beta, s, sigma, gamma, eta0, 0.1, 100, verb); +} + + +struct randwellgetgrad { + randwellgetgrad(double c, unsigned n, unsigned ord, const gsl_vector *k, const gsl_vector *a, const gsl_vector *phi, double w, double s): c(c), n(n), ord(ord), k(k), a(a), phi(phi), w(w), s(s) {} + void operator()(gsl_vector* grad, unsigned size, gsl_vector *state) {domain_energy_randWellGradient(grad, n, state, c, ord, k, a, phi, w, s);} + + private: + double c; + double s; + double w; + unsigned n; + unsigned ord; + const gsl_vector *k; + const gsl_vector *a; + const gsl_vector *phi; +}; + +struct randwellgethess { + randwellgethess(double c, unsigned n, unsigned ord, const gsl_vector *k, const gsl_vector *a, const gsl_vector *phi, double w, double s): c(c), n(n), ord(ord), k(k), a(a), phi(phi), w(w), s(s) {} + void operator()(gsl_matrix* hess, unsigned size, gsl_vector *state) {domain_energy_randWellHessian(hess, n, state, c, ord, k, a, phi, w, s);} + + private: + double c; + double s; + double w; + unsigned n; + unsigned ord; + const gsl_vector *k; + const gsl_vector *a; + const gsl_vector *phi; +}; + +struct randwellgetenergy { + randwellgetenergy(double c, unsigned n, unsigned ord, const gsl_vector *k, const gsl_vector *a, const gsl_vector *phi, double w, double s): c(c), n(n), ord(ord), k(k), a(a), phi(phi), w(w), s(s) {} + double operator()(unsigned size, gsl_vector *state) {return domain_energy_randWellLagrangian(n, state, c, ord, k, a, phi, w, s);} + + private: + double c; + double s; + double w; + unsigned n; + unsigned ord; + const gsl_vector *k; + const gsl_vector *a; + const gsl_vector *phi; +}; + +// Carries out Newton's method. +int domain_minimize_randWell(gsl_vector *z, unsigned n, double c, unsigned ord, const gsl_vector *k, const gsl_vector *a, const gsl_vector *phi, double w, double ss, double eps, unsigned N, double beta, double s, double sigma, double gamma, double eta0, bool verb) { + + unsigned size = 3 * n + 2; + unsigned params = 2 * n + 1; + randwellgetgrad grad(c, n, ord, k, a, phi, w, ss); + randwellgethess hess(c, n, ord, k, a, phi, w, ss); + randwellgetenergy energy(c, n, ord, k, a, phi, w, ss); + + return domain_newton(z, size, params, energy, grad, hess, eps, N, beta, s, sigma, gamma, eta0, 0.1, 100, verb); +} + + +struct fixedmingetgrad { + fixedmingetgrad(double c, unsigned n): c(c), n(n) {} + void operator()(gsl_vector* grad, unsigned size, gsl_vector *state) {domain_energy_fixedGradient(grad, n, state, c);} + + private: + double c; + unsigned n; +}; + +struct fixedmingethess { + fixedmingethess(double c, unsigned n): c(c), n(n) {} + void operator()(gsl_matrix* hess, unsigned size, gsl_vector *state) {domain_energy_fixedHessian(hess, n, state, c);} + + private: + double c; + unsigned n; +}; + +struct fixedmingetenergy { + fixedmingetenergy(double c, unsigned n): c(c), n(n) {} + double operator()(unsigned size, gsl_vector *state) {return domain_energy_fixedLagrangian(n, state, c);} + + private: + double c; + unsigned n; +}; + +// Carries out Newton's method. +int domain_minimize_fixedmin(gsl_vector *z, unsigned n, double c, double eps, unsigned N, double beta, double s, double sigma, double gamma, double bound, bool verb) { + + unsigned size = 3 * n + 3; + unsigned params = 2 * n; + fixedmingetgrad grad(c, n); + fixedmingethess hess(c, n); + fixedmingetenergy energy(c, n); + + return domain_newton(z, size, params, energy, grad, hess, eps, N, beta, s, sigma, gamma, bound, 0.1, 10, verb); +} diff --git a/src/domain_minimize.h b/src/domain_minimize.h new file mode 100644 index 0000000..3754e6a --- /dev/null +++ b/src/domain_minimize.h @@ -0,0 +1,24 @@ +#ifndef DOMAIN_MINIMIZE_H +#define DOMAIN_MINIMIZE_H + +// GSL includes. +#include <gsl/gsl_vector.h> + +// Gives the necessary functions for the Lagrangian, gradient, and Hessian. +#include "domain_energy.h" +#include "domain_newton.h" + +int domain_minimize_naked(gsl_vector *z, unsigned n, double c, double eps, unsigned N, double beta, double s, double sigma, double gamma, double eta0, bool verb); + +int domain_minimize_fixed(gsl_vector *z, unsigned n, double c, double eps, unsigned N, double beta, double s, double sigma); + +int domain_minimize_rand(gsl_vector *z, unsigned n, double c, unsigned ord, const gsl_vector *k, const gsl_vector *a, const gsl_vector *phi, double eps, unsigned N, double beta, double s, double sigma, double gamma, double bound, bool verb); + +int domain_minimize_nakedWell(gsl_vector *z, unsigned n, double c, double w, double ss, double eps, unsigned N, double beta, double s, double sigma, double gamma, double bound, bool verb); + +int domain_minimize_randWell(gsl_vector *z, unsigned n, double c, unsigned ord, const gsl_vector *k, const gsl_vector *a, const gsl_vector *phi, double w, double ss, double eps, unsigned N, double beta, double s, double sigma, double gamma, double bound, bool verb); + +int domain_minimize_fixedmin(gsl_vector *z, unsigned n, double c, double eps, unsigned N, double beta, double s, double sigma, double gamma, double bound, bool verb); + +#endif + diff --git a/src/domain_newton.h b/src/domain_newton.h new file mode 100644 index 0000000..abb7ea3 --- /dev/null +++ b/src/domain_newton.h @@ -0,0 +1,224 @@ +#ifndef DOMAIN_NEWTON_H +#define DOMAIN_NEWTON_H + +#include <math.h> + +// GSL includes. +#include <gsl/gsl_math.h> +#include <gsl/gsl_vector.h> +#include <gsl/gsl_blas.h> +#include <gsl/gsl_sf.h> + +// Eigen's linear solving uses cheap parallelization. +#include <eigen3/Eigen/Dense> + +/* This function is templated so that any set of functions which return an + * energy, gradient, and Hessian given an empty object, the size of the state + * vector, and the state vector can be used. This allows many such sets of + * functions, e.g., that for a fixed domain or a domain on a random background, + * to be used. See the file domain_minimize.cpp for examples of construction + * of these functions. + */ +template <class energy_func, class grad_func, class hess_func> + +int domain_newton(gsl_vector *state, unsigned size, unsigned params, + energy_func get_energy, grad_func get_grad, hess_func get_hess, double + epsilon, unsigned max_iterations, double beta, double s, double sigma, + double gamma, double eta_0, double delta, double bound, bool verbose) { +/* The function domain_newton carries out a modified version of Newton's + * method. On success, 0 is returned. On failure, 1 is returned. + * + * state - GSL_VECTOR + * On entry, state gives the system's initial condition. On + * exit, state contains the result Newton's method. + * + * size - UNSIGNED INTEGER + * On entry, size gives the size of the vector state. Unchanged + * on exit. + * + * params - UNSIGNED INTEGER + * On entry, params gives the number of non-multiplier elements + * in state, which are assumed by the function to be the first + * elements of state. Unchanged on exit. + * + * get_energy - ENERGY_FUNC + * On entry, get_energy is a function that returns a double + * float. The first argument of get_energy is an unsigned + * integer and the second argument is a gsl_vector object. This + * function is expected to take size and state, respectively, + * and return the energy of that state. Unchanged on exit. + * + * get_grad - GRAD_FUNC + * On entry, get_grad is a function that returns void. The + * first argument of get_grad is a gsl_vector object, the second + * argument of get_grad is an unsigned integer, and the third + * argument of get_grad is a gsl_vector object. This function + * is expected to take a vector of size size, size, and state, + * respectively. It leaves the gradient of the energy function + * in the first argument. Unchanged on exit. + * + * get_hess - HESS_FUNC + * On entry, get_hess is a function that returns void. The + * first argument of get_hess is a gsl_matrix object, the second + * argument of get_hess is an unsigned integer, and the third + * argument of get_hess is a gsl_vector object. This function + * is expected to take a matrix of size size by size, size, and + * state, respectively. It leaves the Hessian of the energy + * function in the first argument. Unchanged on exit. + * + * epsilon - DOUBLE FLOAT + * On entry, epsilon gives the number that is used to judge + * convergence. When the norm of the gradient is less than + * epsilon * size, the process is deemed complete and the + * iterations are stopped. Unchanged on exit. + * + * max_iterations - UNSIGNED INTEGER + * On entry, max_iterations gives the maximum number of times + * the algorithm will repeat before failing. Unchanged on exit. + * + * beta - DOUBLE FLOAT + * On entry, beta gives the number which is exponentiated to + * scale the step size in Newton's method. Unchanged on exit. + * + * s - DOUBLE FLOAT + * On entry, s gives a constant scaling of the step size in + * Newton's method. Unchanged on exit. + * + * sigma - DOUBLE FLOAT + * On entry, sigma gives a scaling to the condition on the step + * size in Newton's method. Unchanged on exit. + * + * gamma - DOUBLE FLOAT + * On entry, gamma gives the amount by which the norm of the + * gradient must change for eta to decrement by a factor delta. + * Unchanged on exit. + * + * eta_0 - DOUBLE FLOAT + * On entry, eta_0 gives the starting value of eta. Unchanged + * on exit. + * + * delta - DOUBLE FLOAT + * On entry, delta gives the factor by which eta is decremented. + * Unchanged on exit. + * + * bound - DOUBLE FLOAT + * On entry, delta gives an upper bound to the gradient norm. + * If surpassed, the execution is halted and the program returns + * failure. Unchanged on exit. + * + * verbose - BOOLEAN + * On entry, verbose indicates whether verbose output will be + * printed to stdout by this program. Unchanged on exit. + */ + + // Declaring variables. + double ratio, norm, old_norm, old_energy, energy, grad_dz_prod, alpha, eta; + unsigned iterations, m; + bool converged, bound_exceeded; + + // Declaring GSL variables. + gsl_vector *grad, *dz; + gsl_matrix *hess; + + // Allocating memory for GSL objects + grad = gsl_vector_alloc(size); + dz = gsl_vector_alloc(size); + hess = gsl_matrix_alloc(size, size); + + // Declaring Eigen map objects to wrap the GSL ones. + Eigen::Map<Eigen::VectorXd> grad_eigen(grad->data, size); + Eigen::Map<Eigen::VectorXd> dz_eigen(dz->data, size); + Eigen::Map<Eigen::MatrixXd> hess_eigen(hess->data, size, size); + + // If epsilon > 0, use its value. Otherwise, set to machine precision. + if (epsilon == 0) epsilon = DBL_EPSILON; + + // Initializes the starting value of old_norm at effectively infinity. + old_norm = 1 / DBL_EPSILON; + + // Start the loop parameter at zero. + iterations = 0; + + /* If the loop ends and this boolean has not been flipped, the program will + * know it has not converged. + */ + converged = false; + + // Initializes the value of eta. + eta = eta_0; + + // Begins the algorithm's loop. + while (iterations < max_iterations) { + + // Gets the energy, gradient and Hessian for this iteration. + old_energy = get_energy(size, state); + get_grad(grad, size, state); + get_hess(hess, size, state); + + // Adds eta along the diagonal of the Hessian for non-multiplier entries. + for (unsigned i = 0; i < params; i++) { + gsl_matrix_set(hess, i, i, gsl_matrix_get(hess, i, i) + eta); + } + + // Use LU decomposition to solve for the next step in Newton's method. + dz_eigen = hess_eigen.lu().solve(grad_eigen); + + // Dots the gradient into the step in order to judge the step size. + gsl_blas_ddot(grad, dz, &grad_dz_prod); + + // Initializes the Armijo counter. + m = 0; + + // This loop determines the Armijo step size. + while (true) { + alpha = gsl_sf_pow_int(beta, m) * s; + gsl_vector_scale(dz, alpha); + gsl_vector_sub(state, dz); + + energy = get_energy(size, state); + + if (fabs(old_energy - energy) >= sigma * alpha * grad_dz_prod) break; + else { + gsl_vector_add(state, dz); + gsl_vector_scale(dz, 1 / alpha); + m++; + } + } + + // Gets the new norm of the gradient for comparison. + norm = gsl_blas_dnrm2(grad) / size; + + // Judges if the norm has changed sufficiently little to decrement eta. + if (fabs(norm - old_norm) < gamma * eta) eta *= delta; + + // Prints several useful statistics for debugging purposes. + if (verbose) printf("m was %i, grad_norm was %e, eta was %e, energy was %e\n", + m, norm, eta, energy); + + // Determines if the process has converged to acceptable precision. + if (norm < epsilon) { + converged = true; + break; + } + + // Causes the program to fail if norm has diverged to a large number. + if (norm > bound) break; + + // Reset the norm for the next iteration. + old_norm = norm; + + // Increment the counter. + iterations++; + } + + // Gotta live free, die hard. No one likes memory leaks. + gsl_vector_free(grad); + gsl_vector_free(dz); + gsl_matrix_free(hess); + + // Return conditions to indicate success or failure. + if (converged) return 0; + else return 1; +} + +#endif diff --git a/src/eigenvalues.cpp b/src/eigenvalues.cpp new file mode 100644 index 0000000..c812d0f --- /dev/null +++ b/src/eigenvalues.cpp @@ -0,0 +1,156 @@ +/* eigenvalues.cpp + * + * Copyright (C) 2013 Jaron Kent-Dobias + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +/* This program allows for the generalized eigenvalues of a modulated domain to + * be computed and returned. + */ + +#include "domain_energy.h" +#include "domain_eigen.h" + +#include <unistd.h> +#include <stdio.h> +#include <iostream> +#include <stdlib.h> +#include <math.h> +#include <string> + +// GSL includes. +#include <gsl/gsl_sf.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_complex.h> +#include <gsl/gsl_complex_math.h> +#include <gsl/gsl_vector.h> +#include <gsl/gsl_permutation.h> +#include <gsl/gsl_permute_vector.h> +#include <gsl/gsl_blas.h> +#include <gsl/gsl_linalg.h> +#include <gsl/gsl_sort_vector.h> + + + +// Initializes the program. +int main(int argc, char *argv[]) { + + int opt, min_fails, eigen_follow, eigen_num, examining; + unsigned n, N, ord, size, params, j, M; + double d, c, dc0, dc, g0, g, eigen_thres, approach_thres, eps, eps2, state, old_state, h, bound; + char *in_filename, *out_filename, *k_filename, *a_filename, *phi_filename, str[19], in; + bool subcrit, reset, rand, verbose, fixed; + + // Setting default values. + + gsl_vector *z, *k, *a, *phi, *old_z, *eigenvalues; + gsl_permutation *eigenorder; + gsl_matrix *hess; + rand = false; + fixed = false; + verbose = false; + N=25; + + while ((opt = getopt(argc, argv, "n:c:i:o:O:K:A:P:e:g:N:b:rvd:M:f")) != -1) { + switch (opt) { + case 'n': + n = atoi(optarg); + break; + case 'c': + c = atof(optarg); + break; + case 'i': + in_filename = optarg; + break; + case 'O': + ord = atoi(optarg); + break; + case 'K': + k_filename = optarg; + break; + case 'A': + a_filename = optarg; + break; + case 'P': + phi_filename = optarg; + break; + case 'N': + N = atoi(optarg); + break; + case 'r': + rand = true; + break; + case 'f': + fixed = true; + break; + default: + exit(EXIT_FAILURE); + } + } + + if (rand || !fixed) { + size = 3 * n + 2; + params = 2 * n + 1; + } else { + size = 3 * n + 3; + params = 2 * n; + } + + z = gsl_vector_alloc(size); + eigenvalues = gsl_vector_alloc(size); + eigenorder = gsl_permutation_alloc(size); + hess = gsl_matrix_alloc(size, size); + if (rand) { + k = gsl_vector_alloc(2 * ord); + a = gsl_vector_alloc(ord); + phi = gsl_vector_alloc(ord); + } + + FILE *in_file = fopen(in_filename, "r"); + gsl_vector_fscanf(in_file, z); + fclose(in_file); + + if (rand) { + FILE *k_file = fopen(k_filename, "r"); + gsl_vector_fscanf(k_file, k); + fclose(k_file); + + FILE *a_file = fopen(a_filename, "r"); + gsl_vector_fscanf(a_file, a); + fclose(a_file); + + FILE *phi_file = fopen(phi_filename, "r"); + gsl_vector_fscanf(phi_file, phi); + fclose(phi_file); + } + + if (rand) domain_energy_nakedRandHessian(hess, n, z, c, ord, k, a, phi); + else { + if (fixed) domain_energy_fixedHessian(hess, n, z, c); + else domain_energy_nakedHessian(hess, n, z, c); + } + + domain_eigen_values(eigenvalues, size, params, hess); + domain_eigen_sort(eigenorder, size, N, eigenvalues); + + for (unsigned i = 0; i < N; i++) { + printf("%e\t", gsl_vector_get(eigenvalues, gsl_permutation_get(eigenorder, i))); + } + printf("\n"); + + gsl_vector_free(z); + + +} diff --git a/src/evolve.cpp b/src/evolve.cpp new file mode 100644 index 0000000..f23e15e --- /dev/null +++ b/src/evolve.cpp @@ -0,0 +1,237 @@ +/* evolve.cpp + * + * Copyright (C) 2013 Jaron Kent-Dobias + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +/* This program allows for the iterated minimization of modulated domains as + * the dimensionless parameter Lambda is varied. + */ + +#include "domain_energy.h" +#include "domain_minimize.h" + +#include <unistd.h> +#include <stdio.h> +#include <iostream> +#include <stdlib.h> +#include <math.h> +#include <string> + +// GSL includes. +#include <gsl/gsl_sf.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_complex.h> +#include <gsl/gsl_complex_math.h> +#include <gsl/gsl_vector.h> +#include <gsl/gsl_permutation.h> +#include <gsl/gsl_permute_vector.h> +#include <gsl/gsl_blas.h> +#include <gsl/gsl_linalg.h> +#include <gsl/gsl_sort_vector.h> + + + +// Initializes the program. +int main(int argc, char *argv[]) { + + int opt, min_fails, eigen_follow, eigen_num, examining; + unsigned n, N, ord, size, params, j, M; + double d, c, dc0, dc, g0, g, eigen_thres, approach_thres, eps, eps2, state, old_state, h, bound, da, w, ss; + char *in_filename, *out_filename, *k_filename, *a_filename, *phi_filename, str[19], in; + bool subcrit, reset, rand, verbose, fixed, well; + + // Setting default values. + + gsl_vector *z, *k, *a, *phi, *old_z; + rand = false; + fixed = false; + well = false; + verbose = false; + j=0; + ss=1; + + while ((opt = getopt(argc, argv, "n:c:i:o:O:K:A:P:e:g:N:b:rvd:M:a:fws:W:j:")) != -1) { + switch (opt) { + case 'n': + n = atoi(optarg); + break; + case 'b': + bound = atof(optarg); + break; + case 'c': + c = atof(optarg); + break; + case 'i': + in_filename = optarg; + break; + case 'o': + out_filename = optarg; + break; + case 'O': + ord = atoi(optarg); + break; + case 'K': + k_filename = optarg; + break; + case 'A': + a_filename = optarg; + break; + case 'P': + phi_filename = optarg; + break; + case 'g': + g0 = atof(optarg); + break; + case 'N': + N = atoi(optarg); + break; + case 'j': + j = atoi(optarg); + break; + case 'M': + M = atoi(optarg); + break; + case 'e': + eps = atof(optarg); + break; + case 'd': + dc0 = atof(optarg); + break; + case 'a': + da = atof(optarg); + break; + case 'r': + rand = true; + break; + case 'f': + fixed = true; + break; + case 'w': + well = true; + break; + case 'W': + w = atof(optarg); + break; + case 's': + ss = atof(optarg); + break; + case 'v': + verbose = true; + break; + default: + exit(EXIT_FAILURE); + } + } + + if (rand || !fixed) { + size = 3 * n + 2; + params = 2 * n + 1; + } else { + size = 3 * n + 3; + params = 2 * n; + } + + z = gsl_vector_alloc(size); + old_z = gsl_vector_alloc(size); + if (rand) { + k = gsl_vector_alloc(2 * ord); + a = gsl_vector_alloc(ord); + phi = gsl_vector_alloc(ord); + } + + FILE *in_file = fopen(in_filename, "r"); + gsl_vector_fscanf(in_file, z); + fclose(in_file); + + if (rand) { + FILE *k_file = fopen(k_filename, "r"); + gsl_vector_fscanf(k_file, k); + fclose(k_file); + + FILE *a_file = fopen(a_filename, "r"); + gsl_vector_fscanf(a_file, a); + fclose(a_file); + + FILE *phi_file = fopen(phi_filename, "r"); + gsl_vector_fscanf(phi_file, phi); + fclose(phi_file); + } + + g = g0; + dc = dc0; + + double beta = 0.9; + double s = 1; + double sigma = 0.5; + + if (rand && well) min_fails = domain_minimize_randWell(z, n, c, ord, k, a, phi, w, ss, eps, N, beta, s, sigma, g, bound, verbose); + else if (rand) min_fails = domain_minimize_rand(z, n, c, ord, k, a, phi, eps, N, beta, s, sigma, g, bound, verbose); + else { + if (fixed) min_fails = domain_minimize_fixedmin(z, n, c, eps, N, beta, ss, sigma, g, bound, verbose); + else { + if (well) min_fails = domain_minimize_nakedWell(z, n, c, w, ss, eps, N, beta, s, sigma, g, bound, verbose); + else min_fails = domain_minimize_naked(z, n, c, eps, N, beta, s, sigma, g, bound, verbose); + } + } + + if (min_fails) { + printf("BIFUR: Initial relaxation failed, exiting.\n"); + FILE *out_file = fopen(out_filename, "w"); + gsl_vector_fprintf(out_file, z, "%.10e"); + fclose(out_file); + return 1; + } + + + while (j < M) { + j += 1; + c += dc; + g = g0; + if (rand) gsl_vector_scale(a, da); + + gsl_vector_memcpy(old_z, z); + + printf("EVOLVE: Step %05d, starting with c = %f.\n", j, c); + + while (true) { + if (rand && well) min_fails = domain_minimize_randWell(z, n, c, ord, k, a, phi, w, ss, eps, N, beta, s, sigma, g, bound, verbose); + else if (rand) min_fails = domain_minimize_rand(z, n, c, ord, k, a, phi, eps, N, beta, s, sigma, g, bound, verbose); + else if (fixed) min_fails = domain_minimize_fixedmin(z, n, c, eps, N, beta, s, sigma, g, bound, verbose); + else if (well) min_fails = domain_minimize_nakedWell(z, n, c, w, ss, eps, N, beta, s, sigma, g, bound, verbose); + else min_fails = domain_minimize_naked(z, n, c, eps, N, beta, s, sigma, g, bound, verbose); + + if (!min_fails) break; + printf("EVOLVE: Newton's method failed to converge, reducing gamma.\n"); + gsl_vector_memcpy(z, old_z); + g *= 0.1; + } + + sprintf(str, "output/out-%05d.dat", j); + FILE *fout = fopen(str, "w"); + fprintf(fout, "%.10e\n", c); + gsl_vector_fprintf(fout, z, "%.10e"); + fclose(fout); + } + + FILE *out_file = fopen(out_filename, "w"); + gsl_vector_fprintf(out_file, z, "%.10e"); + fclose(out_file); + + gsl_vector_free(z); + + return 0; + +} diff --git a/src/gradget.cpp b/src/gradget.cpp new file mode 100644 index 0000000..2f8128f --- /dev/null +++ b/src/gradget.cpp @@ -0,0 +1,107 @@ +/* eigenget.cpp + * + * Copyright (C) 2013 Jaron Kent-Dobias + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +// A utility which returns a given generalized eigenvector to the user. + +#include "domain_energy.h" +#include "domain_eigen.h" + +#include <unistd.h> +#include <stdio.h> +#include <iostream> +#include <stdlib.h> +#include <math.h> +#include <string> + +// GSL includes. +#include <gsl/gsl_sf.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_complex.h> +#include <gsl/gsl_complex_math.h> +#include <gsl/gsl_vector.h> +#include <gsl/gsl_permutation.h> +#include <gsl/gsl_permute_vector.h> +#include <gsl/gsl_blas.h> +#include <gsl/gsl_linalg.h> +#include <gsl/gsl_sort_vector.h> + + +// Initializes the program. +int main(int argc, char *argv[]) { + + int opt; + unsigned n, m, nn; + double c, d; + char *in_filename, *out_filename; + bool truehessian; + + gsl_vector *z, *grad, *testgrad; + + truehessian = false; + + while ((opt = getopt(argc, argv, "n:m:c:i:o:td:")) != -1) { + switch (opt) { + case 'n': + n = atoi(optarg); + break; + case 'm': + m = atoi(optarg); + break; + case 'c': + c = atof(optarg); + break; + case 'd': + d = atof(optarg); + break; + case 'i': + in_filename = optarg; + break; + case 'o': + out_filename = optarg; + break; + case 't': + truehessian = true; + break; + default: + exit(EXIT_FAILURE); + } + } + + if (truehessian) nn = 3 * n + 3; + else nn = 3 * n + 3; + + z = gsl_vector_alloc(3 * n + 3); + grad = gsl_vector_alloc(nn); + testgrad = gsl_vector_alloc(nn); + + FILE *in_file = fopen(in_filename, "r"); + gsl_vector_fscanf(in_file, z); + fclose(in_file); + + if (truehessian) domain_energy_fixedGradient(grad, n, z, c); + else domain_energy_fixedGradient(grad, n, z, c); + domain_energy_fixedGradient(testgrad, n, z, c); + + FILE *out_file = fopen(out_filename, "w"); + for (unsigned i = 0; i < nn; i++) { + fprintf(out_file, "%e\t", gsl_vector_get(grad, i) - gsl_vector_get(testgrad, i)); + } + fclose(out_file); +} + + diff --git a/src/hessget.cpp b/src/hessget.cpp new file mode 100644 index 0000000..ad4ae2d --- /dev/null +++ b/src/hessget.cpp @@ -0,0 +1,107 @@ +/* eigenget.cpp + * + * Copyright (C) 2013 Jaron Kent-Dobias + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +// A utility which returns a given generalized eigenvector to the user. + +#include "domain_energy.h" +#include "domain_eigen.h" + +#include <unistd.h> +#include <stdio.h> +#include <iostream> +#include <stdlib.h> +#include <math.h> +#include <string> + +// GSL includes. +#include <gsl/gsl_sf.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_complex.h> +#include <gsl/gsl_complex_math.h> +#include <gsl/gsl_vector.h> +#include <gsl/gsl_permutation.h> +#include <gsl/gsl_permute_vector.h> +#include <gsl/gsl_blas.h> +#include <gsl/gsl_linalg.h> +#include <gsl/gsl_sort_vector.h> + + +// Initializes the program. +int main(int argc, char *argv[]) { + + int opt; + unsigned n, m, nn; + double c; + char *in_filename, *out_filename; + bool truehessian; + + gsl_vector *z; + gsl_matrix *hess, *testhess; + + truehessian = false; + + while ((opt = getopt(argc, argv, "n:m:c:i:o:t")) != -1) { + switch (opt) { + case 'n': + n = atoi(optarg); + break; + case 'm': + m = atoi(optarg); + break; + case 'c': + c = atof(optarg); + break; + case 'i': + in_filename = optarg; + break; + case 'o': + out_filename = optarg; + break; + case 't': + truehessian = true; + break; + default: + exit(EXIT_FAILURE); + } + } + + if (truehessian) nn = 3 * n + 4; + else nn = 3 * n + 3; + + z = gsl_vector_alloc(3 * n + 3); + hess = gsl_matrix_alloc(nn, nn); + testhess = gsl_matrix_alloc(nn, nn); + + FILE *in_file = fopen(in_filename, "r"); + gsl_vector_fscanf(in_file, z); + fclose(in_file); + + domain_energy_fixedHessian(hess, n, z, c); + domain_energy_fixedHessian(testhess, n, z, c); + + FILE *out_file = fopen(out_filename, "w"); + for (unsigned i = 0; i < nn; i++) { + for (unsigned j = 0; j < nn; j++) { + fprintf(out_file, "%e\t", gsl_matrix_get(hess, i, j) - gsl_matrix_get(testhess, i, j)); + } + fprintf(out_file, "\n"); + } + fclose(out_file); +} + + diff --git a/src/initialize.cpp b/src/initialize.cpp new file mode 100644 index 0000000..3c77fc2 --- /dev/null +++ b/src/initialize.cpp @@ -0,0 +1,201 @@ +/* initialize.cpp + * + * Copyright (C) 2013 Jaron Kent-Dobias + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +// Initializes modulated domains. + +#include <unistd.h> +#include <iostream> +#include <string> +#include <stdlib.h> +#include <math.h> +#include <time.h> + +#include <gsl/gsl_sf.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_rng.h> +#include <gsl/gsl_vector.h> + + +int main(int argc, char *argv[]) { + + // Declaring variables. + gsl_vector *z, *k, *a, *phi; + int opt; + unsigned n, m, size, params, word, ord; + double p, k0, a0, R, w0, slope; + char *out_filename, *k_filename, *a_filename, *phi_filename; + bool rand; + + // Default values. + n = 100; + p = 0; + m = 0; + rand = false; + k0 = 1; + a0 = 1; + ord = 0; + word = 0; + w0=1; + slope=1; + + // GNU getopt in action. + while ((opt = getopt(argc, argv, "n:p:m:o:O:rK:A:P:k:a:wW:u:s:")) != -1) { + switch (opt) { + case 'n': + n = atoi(optarg); + break; + case 'p': + p = atof(optarg); + break; + case 'a': + a0 = atof(optarg); + break; + case 'k': + k0 = atof(optarg); + break; + case 'm': + m = atoi(optarg); + break; + case 'O': + ord = atoi(optarg); + break; + case 'o': + out_filename = optarg; + break; + case 'K': + k_filename = optarg; + break; + case 'A': + a_filename = optarg; + break; + case 'P': + phi_filename = optarg; + break; + case 'r': + rand = true; + break; + case 'W': + word = atoi(optarg); + break; + case 'u': + w0 = atof(optarg); + break; + case 's': + slope = atof(optarg); + break; + default: + exit(EXIT_FAILURE); + } + } + + if (rand) { + size = 3 * n + 2; + params = 2 * n + 1; + } else { + size = 3 * n + 3; + params = 2 * n; + } + + z = gsl_vector_alloc(size); + if (rand) { + k = gsl_vector_alloc(2 * (ord + 2 * word)); + a = gsl_vector_alloc(ord + 2 * word); + phi = gsl_vector_alloc(ord + 2* word); + } + + R = sqrt(2 * M_PI / (n * sin(2 * M_PI / n))); + + // setting x[0..n], y[1..n] + for (unsigned i = 0; i < n; i++) { + gsl_vector_set(z, i, R * cos(2 * M_PI * i / n)); + if (rand) gsl_vector_set(z, n + i, R * sin(2 * M_PI * i / n)); + else if (i != 0) gsl_vector_set(z, n + i - 1, R * sin(2 * M_PI * i / n)); + } + + // setting L + gsl_vector_set(z, params - 1, 2 * R * n * sin(M_PI / n)); + + for (unsigned i = 0; i < size - params; i++) { + gsl_vector_set(z, params + i, 1); + } + + if (p > 0) { + for (unsigned i = 0; i < n; i++) { + gsl_vector_set(z, i, gsl_vector_get(z, i) * + (1 + p * cos(2 * M_PI * m * i / n))); + + if (rand) gsl_vector_set(z, n + i, gsl_vector_get(z, n + i) * + (1 + p * cos(2 * M_PI * m * i / n))); + else if (i != 0) gsl_vector_set(z, n + i - 1, gsl_vector_get(z, n + i - 1) * + (1 + p * cos(2 * M_PI * m * i / n))); + } + } + if (rand) { + // The GSL random number generator. Don't try to make more than one random + // background in the same second, or you'll get identical results. + gsl_rng *rand; + rand = gsl_rng_alloc(gsl_rng_ranlux); + gsl_rng_set(rand, time(NULL)); + + double kx, ky, kk; + + for (unsigned i = 0; i < ord; i++) { + while (true) { + kx = gsl_rng_uniform(rand) * 2 - 1; + ky = gsl_rng_uniform(rand) * 2 - 1; + kk = gsl_pow_2(kx)+gsl_pow_2(ky); + + if (kk < 1) { + gsl_vector_set(k, i, k0 * kx); + gsl_vector_set(k, i + ord + 2 * word, k0 * ky); + break; + } + } + gsl_vector_set(a, i, 2 * a0 / ord * gsl_rng_uniform(rand)); + gsl_vector_set(phi, i, 2 * M_PI * gsl_rng_uniform(rand)); + } + } + + if (rand) { + FILE *k_file = fopen(k_filename, "w"); + gsl_vector_fprintf(k_file, k, "%.10g"); + fclose(k_file); + + FILE *a_file = fopen(a_filename, "w"); + gsl_vector_fprintf(a_file, a, "%.10g"); + fclose(a_file); + + FILE *phi_file = fopen(phi_filename, "w"); + gsl_vector_fprintf(phi_file, phi, "%.10g"); + fclose(phi_file); + } + + FILE *out_file = fopen(out_filename, "w"); + gsl_vector_fprintf(out_file, z, "%.10g"); + fclose(out_file); + + + gsl_vector_free(z); + if (rand) { + gsl_vector_free(phi); + gsl_vector_free(k); + gsl_vector_free(a); + } + +} + diff --git a/src/perturb.cpp b/src/perturb.cpp new file mode 100644 index 0000000..9eaf434 --- /dev/null +++ b/src/perturb.cpp @@ -0,0 +1,113 @@ +/* eigenget.cpp + * + * Copyright (C) 2013 Jaron Kent-Dobias + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +// A utility which returns a given generalized eigenvector to the user. + +#include "domain_energy.h" +#include "domain_min.h" +#include "domain_eigen.h" + +#include <unistd.h> +#include <stdio.h> +#include <iostream> +#include <stdlib.h> +#include <math.h> +#include <string> + +// GSL includes. +#include <gsl/gsl_sf.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_complex.h> +#include <gsl/gsl_complex_math.h> +#include <gsl/gsl_vector.h> +#include <gsl/gsl_permutation.h> +#include <gsl/gsl_permute_vector.h> +#include <gsl/gsl_blas.h> +#include <gsl/gsl_linalg.h> +#include <gsl/gsl_sort_vector.h> + + +// Initializes the program. +int main(int argc, char *argv[]) { + + int opt; + unsigned n, m, nn; + double c; + char *in_filename, *out_filename; + bool truehessian; + + truehessian = false; + + gsl_vector *z, *eigenvalues, *eigenvector; + gsl_permutation *eigenorder; + + while ((opt = getopt(argc, argv, "n:m:c:i:o:t")) != -1) { + switch (opt) { + case 'n': + n = atoi(optarg); + break; + case 'm': + m = atoi(optarg); + break; + case 'c': + c = atof(optarg); + break; + case 'i': + in_filename = optarg; + break; + case 'o': + out_filename = optarg; + break; + case 't': + truehessian = true; + break; + default: + exit(EXIT_FAILURE); + } + } + + if (truehessian) nn = 3 * n + 4; + else nn = 3 * n + 3; + + z = gsl_vector_alloc(3 * n + 3); + eigenvalues = gsl_vector_alloc(nn); + eigenvector = gsl_vector_alloc(nn); + eigenorder = gsl_permutation_alloc(nn); + + FILE *in_file = fopen(in_filename, "r"); + gsl_vector_fscanf(in_file, z); + fclose(in_file); + + if (truehessian) { + domain_eigen_truevalues(eigenvalues, n, z, c); + domain_eigen_truesort(eigenorder, n, m, eigenvalues); + domain_eigen_truevector(eigenvector, gsl_permutation_get(eigenorder, m), n, z, c); + } else { + domain_eigen_values(eigenvalues, n, z, c); + domain_eigen_sort(eigenorder, n, m, eigenvalues); + domain_eigen_vector(eigenvector, gsl_permutation_get(eigenorder, m), n, z, c); + } + + FILE *out_file = fopen(out_filename, "w"); + for (unsigned i = 0; i < nn; i++) { + fprintf(out_file, "%e\t", gsl_vector_get(eigenvector, i)); + } + fclose(out_file); +} + + |