sparse-linear-algebra (empty) → 0.1.0.0

raw patch · 9 files changed

+2810/−0 lines, 9 filesdep +QuickCheckdep +basedep +containerssetup-changed

Dependencies added: QuickCheck, base, containers, criterion, hspec, monad-loops, mtl, mwc-random, primitive, sparse-linear-algebra, transformers

Files

• LICENSE +674/−0
• README.md +57/−0
• Setup.hs +2/−0
• app/Main.hs +8/−0
• sparse-linear-algebra.cabal +66/−0
• src/Lib.hs +11/−0
• src/Math/Linear/Sparse.hs +1842/−0
• src/Math/Linear/Sparse/IntMap.hs +149/−0
• test/Spec.hs +1/−0

@@ -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:++    {project}  Copyright (C) {year}  {fullname}+    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>.

@@ -0,0 +1,57 @@+# sparse-linear-algebra++Sparse linear algebra datastructures and algorithms in Haskell++TravisCI : [![Build Status](https://travis-ci.org/ocramz/sparse-linear-algebra.png)](https://travis-ci.org/ocramz/sparse-linear-algebra)++This library provides common numerical analysis functionality, without requiring any external bindings. It is not optimized for performance (yet), but it serves as an experimental platform for scientific computation in a purely functional setting.++Algorithms :++* Iterative linear solvers++    * Conjugate Gradient Squared (CGS)++    * BiConjugate Gradient Stabilized (BiCGSTAB) (non-Hermitian systems)++* Matrix decompositions++    * QR factorization++* Eigenvalue algorithms++    * QR algorithm++    * Rayleigh quotient iteration++* Utilities : Vector and matrix norms, matrix condition number, Givens rotation, Householder reflection++* Predicates : Matrix orthogonality test (A^T A ~= I)++++----------++This is also an experiment in principled scientific programming :++* set the stage by declaring typeclasses and some useful generic operations (normed linear vector spaces, i.e. finite-dimensional spaces equipped with an inner product that induces a distance function),++* define appropriate data structures, and how they relate to those properties (sparse vectors and matrices, defined internally via `Data.IntMap`, are made instances of the VectorSpace and AdditiveGroup classes respectively). This allows to decouple the algorithms from the actual implementation of the backend,++* implement the algorithms, following 1:1 the textbook [1] +++## License++GPL3, see LICENSE++## Credits++Inspired by++* `linear` : https://hackage.haskell.org/package/linear+* `sparse-lin-alg` : https://github.com/laughedelic/sparse-lin-alg++## References++[1] : Y. Saad, Iterative Methods for Sparse Linear Systems, 2nd ed., 2000

@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain

@@ -0,0 +1,8 @@+module Main where++import Lib (ourAdd)++import Text.Printf (printf)++main :: IO ()+main = printf "2 + 3 = %d\n" (ourAdd 2 3)

@@ -0,0 +1,66 @@+name:                sparse-linear-algebra+version:             0.1.0.0+synopsis:            Sparse linear algebra datastructures and algorithms+description:         Please see README.md+homepage:            https://github.com/ocramz/sparse-linear-algebra+license:             BSD3+license-file:        LICENSE+author:              Marco Zocca+maintainer:          zocca.marco gmail+copyright:           2016 Marco Zocca+category:            Math+build-type:          Simple+extra-source-files:  README.md+cabal-version:       >=1.10+tested-with:         GHC == 8.0.1++library+  default-language:    Haskell2010+  ghc-options:         -Wall+  hs-source-dirs:      src+  exposed-modules:     Lib+                       Math.Linear.Sparse+                       Math.Linear.Sparse.IntMap+  build-depends:       QuickCheck+                     , base >= 4.7 && < 5+                     , containers+                     , hspec+                     , primitive >= 0.6.1.0+                     , transformers >= 0.5.2.0+                     -- , lens+                     , mtl >= 2.2.1+                     , mwc-random+                     , monad-loops++executable sparse-linear-algebra+  default-language:    Haskell2010+  ghc-options:         -threaded -rtsopts -with-rtsopts=-N+  hs-source-dirs:      app+  main-is:             Main.hs+  build-depends:       base+                     , mtl >= 2.2.1+                     , mwc-random+                     , primitive >= 0.6.1.0+                     , sparse-linear-algebra+                     , transformers >= 0.5.2.0++test-suite spec+  default-language:    Haskell2010+  ghc-options:         -Wall+  type:                exitcode-stdio-1.0+  hs-source-dirs:      test+  main-is:             Spec.hs+  build-depends:       base+                     , containers+                     , hspec+                     , mtl >= 2.2.1+                     , mwc-random+                     , primitive >= 0.6.1.0+                     , sparse-linear-algebra+                     , transformers >= 0.5.2.0+                     , criterion+                     -- , QuickCheck++source-repository head+  type:     git+  location: https://github.com/ocramz/sparse-linear-algebra

@@ -0,0 +1,11 @@+-- | A library to do stuff.+module Lib+    (+      ourAdd+    ) where++-- | Add two 'Int' values.+ourAdd :: Int  -- ^ left+       -> Int  -- ^ right+       -> Int  -- ^ sum+ourAdd x y = x + y

@@ -0,0 +1,1842 @@+{-# LANGUAGE FlexibleContexts, TypeFamilies, MultiParamTypeClasses, FlexibleInstances #-}+-- {-# OPTIONS_GHC -O2 -rtsopts -with-rtsopts=-K32m -prof#-}++module Math.Linear.Sparse where++import Math.Linear.Sparse.IntMap +++import Control.Monad.Primitive++import Control.Monad (mapM_, forM_, replicateM)+import Control.Monad.Loops++import Control.Monad.Cont+import Control.Monad.State.Strict+import Control.Monad.Writer+import Control.Monad.Trans++import Control.Monad.Trans.State (runStateT)+import Control.Monad.Trans.Writer (runWriterT)++import qualified Data.IntMap.Strict as IM+-- import Data.Utils.StrictFold (foldlStrict) -- hidden in `containers`++import qualified System.Random.MWC as MWC+import qualified System.Random.MWC.Distributions as MWC++import Data.Monoid+import qualified Data.Foldable as F+import qualified Data.Traversable as T+++import qualified Data.List as L+import Data.Maybe++++-- | ========= CLASSES and common operations++-- | Additive ring +class Functor f => Additive f where+  -- | zero element+  zero :: Num a => f a+  +  -- | componentwise operations+  (^+^) :: Num a => f a -> f a -> f a+  (^-^) :: Num a => f a -> f a -> f a++++-- | negate the values in a functor+negated :: (Num a, Functor f) => f a -> f a+negated = fmap negate++x `minus` y = x ^+^ negated y++++-- | Vector space+class Additive f => VectorSpace f where+  -- | multiplication by a scalar+  (.*) :: Num a => a -> f a -> f a+  ++-- |linear interpolation+lerp :: (VectorSpace f, Num a) => a -> f a -> f a -> f a+lerp a u v = a .* u ^+^ ((1-a) .* v)+++-- | Hilbert space (inner product)+class VectorSpace f => Hilbert f where+  -- | inner product+  dot :: Num a => f a -> f a -> a++class Hilbert f => Normed f where+  norm :: (Floating a, Eq a) => a -> f a -> a+++-- some norms and related results++-- squared norm +normSq :: (Hilbert f, Num a) => f a -> a+normSq v = v `dot` v+++-- L1 norm+norm1 :: (Foldable t, Num a, Functor t) => t a -> a+norm1 v = sum (fmap abs v)++-- Euclidean norm+norm2 :: (Hilbert f, Floating a) => f a -> a+norm2 v = sqrt (normSq v)++-- Lp norm (p > 0)+normP :: (Foldable t, Functor t, Floating a) => a -> t a -> a+normP p v = sum u**(1/p) where+  u = fmap (**p) v++-- infinity-norm+normInfty :: (Foldable t, Ord a) => t a -> a+normInfty = maximum++++-- normalize+normalize :: (Normed f, Floating a, Eq a) => a -> f a -> f a+normalize n v = (1 / norm n v) .* v++++++++-- -- Lp inner product (p > 0)+dotLp :: (Set t, Foldable t, Floating a) => a -> t a -> t a ->  a+dotLp p v1 v2 = sum u**(1/p) where+  f a b = (a*b)**p+  u = liftI2 f v1 v2+++-- reciprocal+reciprocal :: (Functor f, Fractional b) => f b -> f b+reciprocal = fmap recip+++-- scale+scale :: (Num b, Functor f) => b -> f b -> f b+scale n = fmap (* n)++++++++-- | FiniteDim : finite-dimensional objects++class Additive f => FiniteDim f where+  type FDSize f :: *+  dim :: f a -> FDSize f+++instance FiniteDim SpVector where+  type FDSize SpVector = Int+  dim = svDim+++instance FiniteDim SpMatrix where+  type FDSize SpMatrix = (Rows, Cols)+  dim = smDim+++-- unary dimension-checking bracket+withDim :: (FiniteDim f, Show e) =>+     f a+     -> (FDSize f -> f a -> Bool)+     -> (f a -> c)+     -> String+     -> (f a -> e)+     -> c+withDim x p f e ef | p (dim x) x = f x+                   | otherwise = error e' where e' = e ++ show (ef x)++-- binary dimension-checking bracket+withDim2 :: (FiniteDim f, FiniteDim g, Show e) =>+     f a+     -> g b+     -> (FDSize f -> FDSize g -> f a -> g b -> Bool)+     -> (f a -> g b -> c)+     -> String+     -> (f a -> g b -> e)+     -> c+withDim2 x y p f e ef | p (dim x) (dim y) x y = f x y+                      | otherwise = error e' where e' = e ++ show (ef x y)+++++++-- | HasData : accessing inner data (do not export)++class Additive f => HasData f a where+  type HDData f a :: * +  dat :: f a -> HDData f a++instance HasData SpVector a where+  type HDData SpVector a = IM.IntMap a+  dat = svData++instance HasData SpMatrix a where+  type HDData SpMatrix a = IM.IntMap (IM.IntMap a)+  dat = smData+++++++-- | Sparse : sparse datastructures++class (FiniteDim f, HasData f a) => Sparse f a where+  spy :: Fractional b => f a -> b+++instance Sparse SpVector a where+  spy = spySV++instance Sparse SpMatrix a where+  spy = spySM+++++class Functor f => Set f where+  -- |union binary lift+  liftU2 :: (a -> a -> a) -> f a -> f a -> f a++  -- |intersection binary lift+  liftI2 :: (a -> b -> c) -> f a -> f b -> f c  ++++-- class (Set f, Sparse f a) => SparseSet f a++-- instance SparseSet SpVector a where++++-- | =======================================================++instance Set IM.IntMap where+  liftU2 = IM.unionWith+  {-# INLINE liftU2 #-}+  liftI2 = IM.intersectionWith+  {-# INLINE liftI2 #-}++-- | IntMap implementation+instance Additive IM.IntMap where+  zero = IM.empty+  {-# INLINE zero #-}+  (^+^) = liftU2 (+)+  {-# INLINE (^+^) #-}+  x ^-^ y = x ^+^ negated y+  {-# INLINE (^-^) #-}++instance VectorSpace IM.IntMap where+  n .* im = IM.map (* n) im+  +instance Hilbert IM.IntMap where+   a `dot` b = sum $ liftI2 (*) a b+              ++instance Normed IM.IntMap where+  norm p v | p==1 = norm1 v+           | p==2 = norm2 v+           | otherwise = normP p v++++-- | =======================================================++-- | Sparse Vector+data SpVector a = SV { svDim :: Int ,+                       svData :: IM.IntMap a} deriving Eq++dimSV :: SpVector a -> Int+dimSV = svDim++spySV :: Fractional b => SpVector a -> b+spySV s = fromIntegral (IM.size (dat s)) / fromIntegral (svDim s)+++-- internal : projection functions, do not export+imSV :: SpVector a -> IM.IntMap a+imSV = svData+++++-- | instances for SpVector+instance Functor SpVector where+  fmap f (SV n x) = SV n (fmap f x)++instance Set SpVector where  +  liftU2 f2 (SV n1 x1) (SV n2 x2) = SV (max n1 n2) (liftU2 f2 x1 x2)+  liftI2 f2 (SV n1 x1) (SV n2 x2) = SV (max n1 n2) (liftI2 f2 x1 x2)+  +instance Foldable SpVector where+    foldr f d v = F.foldr f d (svData v)++instance Additive SpVector where+  zero = SV 0 IM.empty+  (^+^) = liftU2 (+)+  (^-^) = liftU2 (-)++                      +instance VectorSpace SpVector where+  n .* v = scale n v++instance Hilbert SpVector where+  a `dot` b | dim a == dim b = dot (dat a) (dat b)+            | otherwise =+                     error $ "dot : sizes must coincide, instead we got " +++                           show (dim a, dim b)+++instance Normed SpVector where+  norm p (SV _ v) = norm p v+++++++-- | empty sparse vector (size n, no entries)++zeroSV :: Int -> SpVector a+zeroSV n = SV n IM.empty+++singletonSV :: a -> SpVector a+singletonSV x = SV 1 (IM.singleton 0 x)+++-- | create a sparse vector from an association list while discarding all zero entries+mkSpVector :: (Num a, Eq a) => Int -> IM.IntMap a -> SpVector a+mkSpVector d im = SV d $ IM.filterWithKey (\k v -> v /= 0 && inBounds0 d k) im++-- | ", from logically dense array (consecutive indices)+mkSpVectorD :: (Num a, Eq a) => Int -> [a] -> SpVector a+mkSpVectorD d ll = mkSpVector d (IM.fromList $ denseIxArray (take d ll))++-- ", don't filter zero elements+mkSpVector1 :: Int -> IM.IntMap a -> SpVector a+mkSpVector1 d ll = SV d $ IM.filterWithKey (\ k _ -> inBounds0 d k) ll++mkSpVector1D :: Int -> [a] -> SpVector a+mkSpVector1D d ll = mkSpVector1 d (IM.fromList $ denseIxArray (take d ll))++++-- | DENSE vector of `1`s+onesSV :: Num a => Int -> SpVector a+onesSV d = SV d $ IM.fromList $ denseIxArray $ replicate d 1++-- | DENSE vector of `0`s+zerosSV :: Num a => Int -> SpVector a+zerosSV d = SV d $ IM.fromList $ denseIxArray $ replicate d 0+++++-- insert+insertSpVector :: Int -> a -> SpVector a -> SpVector a+insertSpVector i x (SV d xim)+  | inBounds0 d i = SV d (IM.insert i x xim)+  | otherwise = error "insertSpVector : index out of bounds"+++fromListSV :: Int -> [(Int, a)] -> SpVector a+fromListSV d iix = SV d (IM.fromList (filter (inBounds0 d . fst) iix ))++-- toList+toListSV :: SpVector a -> [(IM.Key, a)]+toListSV sv = IM.toList (imSV sv)++-- to dense list (default = 0)+toDenseListSV :: Num b => SpVector b -> [b]+toDenseListSV (SV d im) = fmap (\i -> IM.findWithDefault 0 i im) [0 .. d-1]++++++++++  +instance Show a => Show (SpVector a) where+  show (SV d x) = "SV (" ++ show d ++ ") "++ show (IM.toList x)+++-- | lookup++lookupDenseSV :: Num a => IM.Key -> SpVector a -> a+lookupDenseSV i (SV _ im) = IM.findWithDefault 0 i im ++findWithDefault0IM :: Num a => IM.Key -> IM.IntMap a -> a+findWithDefault0IM = IM.findWithDefault 0+++++-- | SV manipulation++tailSV :: SpVector a -> SpVector a+tailSV (SV n sv) = SV (n-1) ta where+  ta = IM.mapKeys (\i -> i - 1) $ IM.delete 0 sv+  ++headSV :: Num a => SpVector a -> a+headSV sv = fromMaybe 0 (IM.lookup 0 (imSV sv))++++-- | concatenate SpVector+++concatSV :: SpVector a -> SpVector a -> SpVector a+concatSV (SV n1 s1) (SV n2 s2) = SV (n1+n2) (IM.union s1 s2') where+  s2' = IM.mapKeys (+ n1) s2+++++++++++-- | promote a SV to SM++svToSM :: SpVector a -> SpMatrix a+svToSM (SV n d) = SM (n, 1) $ IM.singleton 0 d++++    ++-- | outer vector product++outerProdSV, (><) :: Num a => SpVector a -> SpVector a -> SpMatrix a+outerProdSV v1 v2 = fromListSM (m, n) ixy where+  m = dim v1+  n = dim v2+  ixy = [(i,j, x * y) | (i,x) <- toListSV v1 , (j, y) <- toListSV v2]++(><) = outerProdSV++++++++++-- | =======================================================+++++data SpMatrix a = SM {smDim :: (Rows, Cols),+                      smData :: IM.IntMap (IM.IntMap a)} deriving Eq++++-- | instances for SpMatrix+instance Show a => Show (SpMatrix a) where+  show sm@(SM _ x) = "SM " ++ sizeStr sm ++ " "++ show (IM.toList x)++instance Functor SpMatrix where+  fmap f (SM d md) = SM d ((fmap . fmap) f md)++instance Set SpMatrix where+  liftU2 f2 (SM n1 x1) (SM n2 x2) = SM (maxTup n1 n2) ((liftU2.liftU2) f2 x1 x2)+  liftI2 f2 (SM n1 x1) (SM n2 x2) = SM (minTup n1 n2) ((liftI2.liftI2) f2 x1 x2)+  +instance Additive SpMatrix where+  zero = SM (0,0) IM.empty+  (^+^) = liftU2 (+)+  (^-^) = liftU2 (-)+++-- | TODO : use semilattice properties instead+maxTup, minTup :: Ord t => (t, t) -> (t, t) -> (t, t)+maxTup (x1,y1) (x2,y2) = (max x1 x2, max y1 y2)+minTup (x1,y1) (x2,y2) = (min x1 x2, min y1 y2)++-- | empty matrix of size d+emptySpMatrix :: (Int, Int) -> SpMatrix a+emptySpMatrix d = SM d IM.empty++++-- multiply matrix by a scalar+matScale :: Num a => a -> SpMatrix a -> SpMatrix a+matScale a = fmap (*a)++-- Frobenius norm (sqrt of trace of M^T M)+normFrobenius :: SpMatrix Double -> Double+normFrobenius m = sqrt $ foldlSM (+) 0 m' where+  m' | nrows m > ncols m = transposeSM m ## m+     | otherwise = m ## transposeSM m +  ++++++-- | ========= MATRIX METADATA++-- type synonyms+type Rows = Int+type Cols = Int++type IxRow = Int+type IxCol = Int++-- -- predicates+-- are the supplied indices within matrix bounds?+validIxSM :: SpMatrix a -> (Int, Int) -> Bool+validIxSM mm = inBounds02 (dim mm)++-- is the matrix square?+isSquareSM :: SpMatrix a -> Bool+isSquareSM m = nrows m == ncols m++-- is the matrix diagonal?+isDiagonalSM :: SpMatrix a -> Bool+isDiagonalSM m = IM.size d == nrows m where+  d = IM.filterWithKey ff (immSM m)+  ff irow row = IM.size row == 1 &&+                IM.size (IM.filterWithKey (\j _ -> j == irow) row) == 1+++++++++-- -- internal projection functions, do not export:+immSM :: SpMatrix t -> IM.IntMap (IM.IntMap t)+immSM (SM _ imm) = imm++dimSM :: SpMatrix t -> (Rows, Cols)+dimSM (SM d _) = d++nelSM :: SpMatrix t -> Int+nelSM (SM (nr,nc) _) = nr*nc++-- | nrows, ncols : size accessors+nrows, ncols :: SpMatrix a -> Int+nrows = fst . dim+ncols = snd . dim+++++data SMInfo = SMInfo { smNz :: Int,+                       smSpy :: Double} deriving (Eq, Show)++infoSM :: SpMatrix a -> SMInfo+infoSM s = SMInfo (nzSM s) (spySM s)++nzSM :: SpMatrix a -> Int+nzSM s = sum $ fmap IM.size (immSM s)++spySM :: Fractional b => SpMatrix a -> b+spySM s = fromIntegral (nzSM s) / fromIntegral (nelSM s)+++-- # NZ in row i++nzRowU :: SpMatrix a -> IM.Key -> Int+nzRowU s i = maybe 0 IM.size (IM.lookup i $ immSM s)++nzRow :: SpMatrix a -> IM.Key -> Int+nzRow s i | inBounds0 (nrows s) i = nzRowU s i+          | otherwise = error "nzRow : index out of bounds"+++++-- | bandwidth bounds (min, max)++bwMinSM :: SpMatrix a -> Int+bwMinSM = fst . bwBoundsSM++bwMaxSM :: SpMatrix a -> Int+bwMaxSM = snd . bwBoundsSM++bwBoundsSM :: SpMatrix a -> (Int, Int)+bwBoundsSM s = -- b+                (snd $ IM.findMin b,+                snd $ IM.findMax b)+  where+  ss = immSM s+  fmi = fst . IM.findMin+  fma = fst . IM.findMax+  b = fmap (\x -> fma x - fmi x + 1:: Int) ss++++++++++-- | ========= SPARSE MATRIX BUILDERS++zeroSM :: Int -> Int -> SpMatrix a+zeroSM m n = SM (m,n) IM.empty +++insertSpMatrix :: IxRow -> IxCol -> a -> SpMatrix a -> SpMatrix a+insertSpMatrix i j x s+  | inBounds02 d (i,j) = SM d $ insertIM2 i j x smd +  | otherwise = error "insertSpMatrix : index out of bounds" where+      smd = immSM s+      d = dim s+++-- | from list (row, col, value)+fromListSM' :: Foldable t => t (IxRow, IxCol, a) -> SpMatrix a -> SpMatrix a+fromListSM' iix sm = foldl ins sm iix where+  ins t (i,j,x) = insertSpMatrix i j x t++fromListSM :: Foldable t => (Int, Int) -> t (IxRow, IxCol, a) -> SpMatrix a+fromListSM (m,n) iix = fromListSM' iix (zeroSM m n)+++fromListDenseSM :: Int -> [a] -> SpMatrix a+fromListDenseSM m ll = fromListSM (m, n) $ denseIxArray2 m ll where+  n = length ll `div` m+  +++-- | to List++-- toDenseListSM : populate missing entries with 0+toDenseListSM :: Num t => SpMatrix t -> [(IxRow, IxCol, t)]+toDenseListSM m =+  [(i, j, m @@ (i, j)) | i <- [0 .. nrows m - 1], j <- [0 .. ncols m- 1]]++++++-- -- create diagonal and identity matrix+mkDiagonal :: Int -> [a] -> SpMatrix a+mkDiagonal n = mkSubDiagonal n 0+++eye :: Num a => Int -> SpMatrix a+eye n = mkDiagonal n (ones n)++ones :: Num a => Int -> [a]+ones n = replicate n 1+++  +++-- super- and sub- diagonal++mkSubDiagonal :: Int -> Int -> [a] -> SpMatrix a+mkSubDiagonal n o xx | abs o < n = if o >= 0+                                   then fz ii jj xx+                                   else fz jj ii xx+                     | otherwise = error "mkSubDiagonal : offset > dimension" where+  ii = [0 .. n-1]+  jj = [abs o .. n - 1]+  fz a b x = fromListSM (n,n) (zip3 a b x)+++-- fromList :: [(Key,a)] -> IntMap a+-- fromList xs+--   = foldlStrict ins empty xs+--   where+--     ins t (k,x)  = insert k x t+++++encode :: (Int, Int) -> (Rows, Cols) -> Int+encode (nr,_) (i,j) = i + (j * nr)++decode :: (Int, Int) -> Int -> (Rows, Cols)+decode (nr, _) ci = (r, c) where (c,r ) = quotRem ci nr++++++++++-- | ========= SUB-MATRICES+++extractSubmatrixSM :: SpMatrix a -> (Int, Int) -> (Int, Int) -> SpMatrix a+extractSubmatrixSM (SM (r, c) im) (i1, i2) (j1, j2)+  | q = SM (m', n') imm'+  | otherwise = error $ "extractSubmatrixSM : invalid indexing " ++ show (i1, i2) ++ ", " ++ show (j1, j2) where+  imm' = mapKeysIM2 (\i -> i - i1) (\j -> j - j1) $  -- rebalance keys+          IM.filter (not . IM.null) $                -- remove all-null rows+          ifilterIM2 ff im                           -- keep `submatrix`+  ff i j _ = i1 <= i &&+             i <= i2 &&+             j1 <= j &&+             j <= j2+  (m', n') = (i2-i1 + 1, j2-j1 + 1)+  q = inBounds0 r i1  &&+      inBounds0 r i2 &&+      inBounds0 c j1  &&+      inBounds0 c j2 &&      +      i2 >= i1++-- extract row / column+extractRowSM :: SpMatrix a -> Int -> SpMatrix a+extractRowSM sm i = extractSubmatrixSM sm (i, i) (0, ncols sm - 1)++extractColSM :: SpMatrix a -> Int -> SpMatrix a+extractColSM sm j = extractSubmatrixSM sm (0, nrows sm - 1) (j, j)++++-- demote (n x 1) or (1 x n) SpMatrix to SpVector+toSV :: SpMatrix a -> SpVector a+toSV (SM (m,n) im) = SV d $ snd . head $ IM.toList im where+  d | m==1 && n==1 = 1+    | m==1 && n>1 = n +    | n==1 && m>1 = m+    | otherwise = error $ "toSV : incompatible dimensions " ++ show (m,n)+++-- extract row or column and place into SpVector+extractCol :: SpMatrix a -> Int -> SpVector a+extractCol m i = toSV $ extractColSM m i++extractRow :: SpMatrix a -> Int -> SpVector a+extractRow m j = toSV $ extractRowSM m j++++++++-- | ========= MATRIX STACKING++vertStackSM, (-=-) :: SpMatrix a -> SpMatrix a -> SpMatrix a+vertStackSM mm1 mm2 = SM (m, n) $ IM.union u1 u2 where+  nro1 = nrows mm1+  m = nro1 + nrows mm2+  n = max (ncols mm1) (ncols mm2)+  u1 = immSM mm1+  u2 = IM.mapKeys (+ nro1) (immSM mm2)++(-=-) = vertStackSM+++horizStackSM, (-||-) :: SpMatrix a -> SpMatrix a -> SpMatrix a+horizStackSM mm1 mm2 = t (t mm1 -=- t mm2) where+  t = transposeSM++(-||-) = horizStackSM++++++++++-- | ========= LOOKUP++lookupSM :: SpMatrix a -> IM.Key -> IM.Key -> Maybe a+lookupSM (SM _ im) i j = IM.lookup i im >>= IM.lookup j++-- | Looks up an element in the matrix (if not found, zero is returned)++lookupWD_SM, (@@) :: Num a => SpMatrix a -> (IM.Key, IM.Key) -> a+lookupWD_SM sm (i,j) =+  fromMaybe 0 (lookupSM sm i j)++lookupWD_IM :: Num a => IM.IntMap (IM.IntMap a) -> (IM.Key, IM.Key) -> a+lookupWD_IM im (i,j) = fromMaybe 0 (IM.lookup i im >>= IM.lookup j)++(@@) = lookupWD_SM+++++-- FIXME : to throw an exception or just ignore the out-of-bound access ?++++++++-- | ========= MISC SpMatrix OPERATIONS++foldlSM :: (a -> b -> b) -> b -> SpMatrix a -> b+foldlSM f n (SM _ m)= foldlIM2 f n m++ifoldlSM :: (IM.Key -> IM.Key -> a -> b -> b) -> b -> SpMatrix a -> b+ifoldlSM f n (SM _ m) = ifoldlIM2' f n m+++++++-- | mapping ++++++-- | folding++-- count sub-diagonal nonzeros+countSubdiagonalNZSM :: SpMatrix a -> Int+countSubdiagonalNZSM (SM _ im) = countSubdiagonalNZ im+++-- | filtering++-- extractDiagonalSM :: (Num a, Eq a) => SpMatrix a -> SpVector a+-- extractDiagonalSM (SM (m,n) im) = mkSpVectorD m $ extractDiagonalIM2 im++-- extract with default 0+extractDiagonalDSM :: Num a => SpMatrix a -> SpVector a+extractDiagonalDSM mm = mkSpVector1D n $ foldr ins [] ll  where+  ll = [0 .. n - 1]+  n = nrows mm+  ins i acc = mm@@(i,i) : acc+  +++++  ++--  filtering the index subset that lies below the diagonal++subdiagIndicesSM :: SpMatrix a -> [(IM.Key, IM.Key)]+subdiagIndicesSM (SM _ im) = subdiagIndices im++++++-- | sparsify : remove 0s (!!!)++sparsifyIM2 :: IM.IntMap (IM.IntMap Double) -> IM.IntMap (IM.IntMap Double)+sparsifyIM2 = ifilterIM2 (\_ _ x -> abs x >= eps)++sparsifySM :: SpMatrix Double -> SpMatrix Double+sparsifySM (SM d im) = SM d $ sparsifyIM2 im++++-- | ROUNDING operations (!!!)+                              +roundZeroOneSM :: SpMatrix Double -> SpMatrix Double+roundZeroOneSM (SM d im) = sparsifySM $ SM d $ mapIM2 roundZeroOne im++++++++  ++++-- | ========= ALGEBRAIC PRIMITIVE OPERATIONS+++-- | transpose+++transposeSM, (#^) :: SpMatrix a -> SpMatrix a+transposeSM (SM (m, n) im) = SM (n, m) (transposeIM2 im)++(#^) = transposeSM++++-- | A^T B+(#^#) :: SpMatrix Double -> SpMatrix Double -> SpMatrix Double+a #^# b = transposeSM a #~# b+++-- | A B^T+(##^) :: SpMatrix Double -> SpMatrix Double -> SpMatrix Double+a ##^ b = a #~# transposeSM b+++++++-- | matrix action on a vector++{- +FIXME : matVec is more general than SpVector's :++\m v -> fmap (`dot` v) m+  :: (Normed f1, Num b, Functor f) => f (f1 b) -> f1 b -> f b+-}++++-- matrix on vector+matVec, (#>) :: Num a => SpMatrix a -> SpVector a -> SpVector a+matVec (SM (nr, nc) mdata) (SV n sv)+  | nc == n = SV nr $ fmap (`dot` sv) mdata+  | otherwise = error $ "matVec : mismatching dimensions " ++ show (nc, n)++(#>) = matVec++-- vector on matrix (FIXME : transposes matrix: more costly than `matVec`)+vecMat, (<#) :: Num a => SpVector a -> SpMatrix a -> SpVector a  +vecMat (SV n sv) (SM (nr, nc) mdata)+  | n == nr = SV nc $ fmap (`dot` sv) (transposeIM2 mdata)+  | otherwise = error $ "vecMat : mismatching dimensions " ++ show (n, nr)++(<#) = vecMat+++++++++++-- matVec' mm vv =+--   withDim2 mm vv (\(nro, nco) nv _ _ -> nco == nv) matVecU "matVec : mismatching dimensions"+--    (\ m v -> unwords [show (dim m), show (dim v)])+++-- asdfm ll = unwords (map (show . dim) ll)+++++-- | matrix-matrix product++-- unsafe matMat+matMatU :: Num a => SpMatrix a -> SpMatrix a -> SpMatrix a+matMatU m1 m2 =+  SM (nrows m1, ncols m2) im where+    im = fmap (\vm1 -> (`dot` vm1) <$> transposeIM2 (immSM m2)) (immSM m1)+++-- matMat, (##) :: Num a => SpMatrix a -> SpMatrix a -> SpMatrix a+-- matMat (SM (nr1,nc1) m1) (SM (nr2,nc2) m2)+--   | nc1 == nr2 = SM (nr1, nc2) $+--       fmap (\vm1 -> fmap (`dot` vm1) (transposeIM2 m2)) m1+--   | otherwise = error "matMat : incompatible matrix sizes"++matMat, (##) :: Num a => SpMatrix a -> SpMatrix a -> SpMatrix a+matMat m1 m2+  | c1 == r2 = matMatU m1 m2+  | otherwise = error $ "matMat : incompatible matrix sizes" ++ show (d1, d2) where+      d1@(r1, c1) = dim m1+      d2@(r2, c2) = dim m2+    ++(##) = matMat++-- matMat m1 m2 =+--   withDim2 m1 m2+--     (\(r1,c1) (r2,c2) _ _ -> c1 == r2)+--     matMatU+--     "matMat : incompatible matrix sizes"+--     (\m1 m2 -> unwords [show (dim m1), show (dim m2)])++++++-- | sparsified matrix-matrix product (prunes all elements `x` for which `abs x <= eps`)+matMatSparsified, (#~#)  :: SpMatrix Double -> SpMatrix Double -> SpMatrix Double+matMatSparsified m1 m2 = sparsifySM $ matMat m1 m2++(#~#) = matMatSparsified+++++++-- | ========= predicates++-- is the matrix orthogonal? i.e. Q^t ## Q == I+isOrthogonalSM :: SpMatrix Double -> Bool+isOrthogonalSM sm@(SM (_,n) _) = rsm == eye n where+  rsm = roundZeroOneSM $ transposeSM sm ## sm++++++++++-- | ========= condition number++-- uses the R matrix from the QR factorization+conditionNumberSM :: SpMatrix Double -> Double+conditionNumberSM m | isInfinite kappa = error "Infinite condition number : rank-deficient system"+                    | otherwise = kappa where+  kappa = lmax / lmin+  (_, r) = qr m+  u = extractDiagonalDSM r  -- FIXME : need to extract with default element 0 +  lmax = abs (maximum u)+  lmin = abs (minimum u)++++++++-- | ========= Householder transformation++hhMat :: Num a => a -> SpVector a -> SpMatrix a+hhMat beta x = eye n ^-^ scale beta (x >< x) where+  n = dim x+++-- a vector `x` uniquely defines an orthogonal plane; the Householder operator reflects any point `v` with respect to this plane:+-- v' = (I - 2 x >< x) v +hhRefl :: SpVector Double -> SpMatrix Double+hhRefl = hhMat 2.0++++++++++++-- | ========= Givens rotation matrix+++hypot :: Floating a => a -> a -> a+hypot x y = abs x * (sqrt (1 + y/x)**2)++sign :: (Ord a, Num a) => a -> a+sign x+  | x > 0 = 1+  | x == 0 = 0+  | otherwise = -1 ++givensCoef :: (Ord a, Floating a) => a -> a -> (a, a, a)+givensCoef a b  -- returns (c, s, r) where r = norm (a, b)+  | b==0 = (sign a, 0, abs a)+  | a==0 = (0, sign b, abs b)+  | abs a > abs b = let t = b/a+                        u = sign a * abs ( sqrt (1 + t**2))+                      in (1/u, - t/u, a*u)+  | otherwise = let t = a/b+                    u = sign b * abs ( sqrt (1 + t**2))+                in (t/u, - 1/u, b*u)+++{-+Givens method, row version: choose other row index i' s.t. i' is :+* below the diagonal+* corresponding element is nonzero++QR.C1 ) To zero out entry A(i, j) we must find row k such that A(k, j) is+non-zero but A has zeros in row k for all columns less than j.+-}++givens :: SpMatrix Double -> Int -> Int -> SpMatrix Double+givens mm i j +  | validIxSM mm (i,j) && isSquareSM mm =+       sparsifySM $ fromListSM' [(i,i,c),(j,j,c),(j,i,-s),(i,j,s)] (eye (nrows mm))+  | otherwise = error "givens : indices out of bounds"      +  where+    (c, s, _) = givensCoef a b+    i' = head $ fromMaybe (error $ "givens: no compatible rows for entry " ++ show (i,j)) (candidateRows (immSM mm) i j)+    a = mm @@ (i', j)+    b = mm @@ (i, j)   -- element to zero out++-- is the `k`th the first nonzero column in the row?+firstNonZeroColumn :: IM.IntMap a -> IM.Key -> Bool+firstNonZeroColumn mm k = isJust (IM.lookup k mm) &&+                          isNothing (IM.lookupLT k mm)++-- returns a set of rows {k} that satisfy QR.C1+candidateRows :: IM.IntMap (IM.IntMap a) -> IM.Key -> IM.Key -> Maybe [IM.Key]+candidateRows mm i j | IM.null u = Nothing+                     | otherwise = Just (IM.keys u) where+  u = IM.filterWithKey (\irow row -> irow /= i &&+                                     firstNonZeroColumn row j) mm++++++-- | ========= QR algorithm++{-+applies Givens rotation iteratively to zero out sub-diagonal elements+-}+++qr :: SpMatrix Double -> (SpMatrix Double, SpMatrix Double)+qr mm = (transposeSM qmatt, rmat)  where+  qmatt = F.foldl' (#~#) ee $ gmats mm -- Q^T = (G_n * G_n-1 ... * G_1)+  rmat = qmatt #~# mm                  -- R = Q^T A+  ee = eye (nrows mm)+      +-- Givens matrices in order [G1, G2, .. , G_N ]+gmats :: SpMatrix Double -> [SpMatrix Double]+gmats mm = gm mm (subdiagIndicesSM mm) where+ gm m ((i,j):is) = let g = givens m i j+                   in g : gm (g #~# m) is+ gm _ [] = []++++++-- -- | QR algorithm, state transformer version+-- gmatST0 (m, (i,j):is) = (m', is) where    -- WRONG, possible access to []+--   g = givens m i j                        +--   m' = g #~# m+-- gmatST0 (m, []) = (eye (nrows m), [])++-- gmatST m = gmatST0 (m, subdiagIndicesSM m)++++++-- | ========= Eigenvalues, using QR+++eigsQR :: Int -> SpMatrix Double -> SpVector Double+eigsQR nitermax m = extractDiagonalDSM $ execState (convergtest eigsStep) m where+  eigsStep m = r #~# q where (q, r) = qr m+  convergtest g = modifyInspectN nitermax f g where+    f [m1, m2] = let dm1 = extractDiagonalDSM m1+                     dm2 = extractDiagonalDSM m2+                 in norm2 (dm1 ^-^ dm2) <= eps+++++++-- | ========= Eigenvalues, using Rayleigh iteration++rayleighStep ::+  SpMatrix Double ->+  (SpVector Double, Double) ->+  (SpVector Double, Double)    -- updated estimate of (eigenvector, eigenvalue)+rayleighStep aa (b, mu) = (b', mu') where+  ii = eye (nrows aa)+  nom = (aa ^-^ (mu `matScale` ii)) <\> b+  b' = normalize 2 nom+  mu' = b' `dot` (aa #> b') / (b' `dot` b')++eigRayleigh :: Int                -- max # iterations+     -> SpMatrix Double           -- matrix+     -> (SpVector Double, Double) -- initial guess of (eigenvector, eigenvalue)+     -> (SpVector Double, Double) -- final estimate of (eigenvector, eigenvalue)+eigRayleigh nitermax m = execState (convergtest (rayleighStep m)) where+  convergtest g = modifyInspectN nitermax f g where+    f [(b1, _), (b2, _)] = norm2 (b2 ^-^ b1) <= eps ++++++-- | ========= Householder vector (G & VL Alg. 5.1.1, function `house`)++-- hhV :: (Ord a, Floating a) => SpVector a -> (SpVector a, a)+hhV :: SpVector Double -> (SpVector Double, Double)+hhV x = (v, beta) where+  n = dim x+  tx = tailSV x+  sigma = tx `dot` tx+  vtemp = singletonSV 1 `concatSV` tx+  (v, beta) | sigma <= eps = (vtemp, 0)+            | otherwise = let mu = sqrt (headSV x**2 + sigma)+                              xh = headSV x+                              vh | xh <= 1 = xh - mu+                                 | otherwise = - sigma / (xh + mu)+                              vnew = (1 / vh) .* insertSpVector 0 vh vtemp     +                          in (vnew, 2 * xh**2 / (sigma + vh**2))++                         +++++++-- | ========= SVD++{- Golub & Van Loan, sec 8.6.2 (p 452 segg.)++SVD of A :++* reduce A to upper bidiagonal form B (Alg. 5.4.2)+* compute SVD of B (implicit-shift QR step, Alg. 8.3.2)++-}+++++++++++++++++++-- | =======================================================++-- | LINEAR SOLVERS : solve A x = b++-- | numerical tolerance for e.g. solution convergence+eps :: Double+eps = 1e-8++-- | residual of candidate solution x0+residual :: Num a => SpMatrix a -> SpVector a -> SpVector a -> SpVector a+residual aa b x0 = b ^-^ (aa #> x0)++converged :: SpMatrix Double -> SpVector Double -> SpVector Double -> Bool+converged aa b x0 = normSq (residual aa b x0) <= eps++++-- | CGS++-- | one step of CGS+cgsStep :: SpMatrix Double -> SpVector Double -> CGS -> CGS+cgsStep aa rhat (CGS x r p u) = CGS xj1 rj1 pj1 uj1+  where+  aap = aa #> p+  alphaj = (r `dot` rhat) / (aap `dot` rhat)+  q = u ^-^ (alphaj .* aap)+  xj1 = x ^+^ (alphaj .* (u ^+^ q))  -- updated solution+  rj1 = r ^-^ (alphaj .* (aa #> (u ^+^ q)))-- updated residual+  betaj = (rj1 `dot` rhat) / (r `dot` rhat)+  uj1 = rj1 ^+^ (betaj .* q)+  pj1 = uj1 ^+^ (betaj .* (q ^+^ (betaj .* p)))++data CGS = CGS { _x :: SpVector Double,+                 _r :: SpVector Double,+                 _p :: SpVector Double,+                 _u :: SpVector Double } deriving Eq++cgs ::+  SpMatrix Double ->+  SpVector Double ->+  SpVector Double ->+  SpVector Double ->+  -- Int ->+  CGS+cgs aa b x0 rhat =+  -- execState (replicateM n (modify (cgsStep aa rhat))) cgsInit where+  execState (untilConverged _x (cgsStep aa rhat)) cgsInit where+  r0 = b ^-^ (aa #> x0)    -- residual of initial guess solution+  p0 = r0+  u0 = r0+  cgsInit = CGS x0 r0 p0 u0+++instance Show CGS where+  show (CGS x r p u) = "x = " ++ show x ++ "\n" +++                                "r = " ++ show r ++ "\n" +++                                "p = " ++ show p ++ "\n" +++                                "u = " ++ show u ++ "\n"++++  ++-- | BiCSSTAB++-- _aa :: SpMatrix Double,    -- matrix+-- _b :: SpVector Double,     -- rhs+-- _r0 :: SpVector Double,    -- initial residual+-- _r0hat :: SpVector Double, -- candidate solution: r0hat `dot` r0 >= 0++-- | one step of BiCGSTAB+bicgstabStep :: SpMatrix Double -> SpVector Double -> BICGSTAB -> BICGSTAB+bicgstabStep aa r0hat (BICGSTAB x r p) = BICGSTAB xj1 rj1 pj1 where+  aap = aa #> p+  alphaj = (r `dot` r0hat) / (aap `dot` r0hat)+  sj = r ^-^ (alphaj .* aap)+  aasj = aa #> sj+  omegaj = (aasj `dot` sj) / (aasj `dot` aasj)+  xj1 = x ^+^ (alphaj .* p) ^+^ (omegaj .* sj)+  rj1 = sj ^-^ (omegaj .* aasj)+  betaj = (rj1 `dot` r0hat)/(r `dot` r0hat) * alphaj / omegaj+  pj1 = rj1 ^+^ (betaj .* (p ^-^ (omegaj .* aap)))++data BICGSTAB = BICGSTAB { _xBicgstab :: SpVector Double,+                           _rBicgstab :: SpVector Double,+                           _pBicgstab :: SpVector Double} deriving Eq+++bicgstab+  :: SpMatrix Double+     -> SpVector Double+     -> SpVector Double+     -> SpVector Double+     -- -> Int+     -> BICGSTAB+bicgstab aa b x0 r0hat =+  -- execState (replicateM n (modify (bicgstabStep aa r0hat))) bicgsInit where+  execState (untilConverged _xBicgstab (bicgstabStep aa r0hat)) bicgsInit where+   r0 = b ^-^ (aa #> x0)    -- residual of initial guess solution+   p0 = r0+   bicgsInit = BICGSTAB x0 r0 p0+   -- q (BICGSTAB xi _ _) = nor++instance Show BICGSTAB where+  show (BICGSTAB x r p) = "x = " ++ show x ++ "\n" +++                                "r = " ++ show r ++ "\n" +++                                "p = " ++ show p ++ "\n"++++++++++++-- | ========= LINEAR SOLVERS INTERFACE++data LinSolveMethod = CGS_ | BICGSTAB_ deriving (Eq, Show) ++-- random starting vector+linSolveM ::+  PrimMonad m =>+    LinSolveMethod -> SpMatrix Double -> SpVector Double -> m (SpVector Double)+linSolveM method aa b = do+  let (m, n) = dim aa+      nb     = dim b+  if n /= nb then error "linSolve : operand dimensions mismatch" else do+    x0 <- randVec nb+    case method of CGS_ -> return $ _xBicgstab (bicgstab aa b x0 x0)+                   BICGSTAB_ -> return $ _x (cgs aa b x0 x0)++-- deterministic starting vector (every component at 0.1) +linSolve ::+  LinSolveMethod -> SpMatrix Double -> SpVector Double -> SpVector Double+linSolve method aa b+  | n /= nb = error "linSolve : operand dimensions mismatch"+  | otherwise = solve aa b where+      solve aa' b' | isDiagonalSM aa = (reciprocal aa') #> b'+                   | otherwise = solveWith aa' b' +      solveWith aa' b' = case method of+                                CGS_ ->  _xBicgstab (bicgstab aa' b' x0 x0)+                                BICGSTAB_ -> _x (cgs aa' b' x0 x0)+      x0 = mkSpVectorD n $ replicate n 0.1 +      (m, n) = dim aa+      nb     = dim b++-- <\> : sets default solver method ++(<\>) :: SpMatrix Double -> SpVector Double -> SpVector Double      +(<\>) = linSolve BICGSTAB_ +  ++++++++-- | TODO : if system is poorly conditioned, is it better to warn the user or just switch solvers (e.g. via the pseudoinverse) ?++-- linSolveQR aa b init f1 stepf+--   | isInfinite k = do+--        tell "linSolveQR : rank-deficient system"+--   | otherwise = do+--        solv aa b init+--     where+--      (q, r) = qr aa+--      k = conditionNumberSM r+--      solv aa b init = execState (untilConverged f1 stepf) init+++++++++++-- | ========= PRETTY PRINTING++-- | Show details and contents of sparse matrix++sizeStr :: SpMatrix a -> String+sizeStr sm =+  unwords ["(",show (nrows sm),"rows,",show (ncols sm),"columns ) ,",show nz,"NZ ( sparsity",show sy,")"] where+  (SMInfo nz sy) = infoSM sm +++showNonZero :: (Show a, Num a, Eq a) => a -> String+showNonZero x  = if x == 0 then " " else show x+++toDenseRow :: Num a => SpMatrix a -> IM.Key -> [a]+toDenseRow (SM (_,ncol) im) irow =+  fmap (\icol -> im `lookupWD_IM` (irow,icol)) [0..ncol-1]++toDenseRowClip :: (Show a, Num a) => SpMatrix a -> IM.Key -> Int -> String+toDenseRowClip sm irow ncomax+  | ncols sm > ncomax = unwords (map show h) ++  " ... " ++ show t+  | otherwise = show dr+     where dr = toDenseRow sm irow+           h = take (ncomax - 2) dr+           t = last dr++newline :: IO ()+newline = putStrLn ""++printDenseSM :: (Show t, Num t) => SpMatrix t -> IO ()+printDenseSM sm = do+  newline+  putStrLn $ sizeStr sm+  newline+  printDenseSM' sm 5 5+  newline+  where    +    printDenseSM' :: (Show t, Num t) => SpMatrix t -> Int -> Int -> IO ()+    printDenseSM' sm'@(SM (nr,_) _) nromax ncomax = mapM_ putStrLn rr_' where+      rr_ = map (\i -> toDenseRowClip sm' i ncomax) [0..nr - 1]+      rr_' | nrows sm > nromax = take (nromax - 2) rr_ ++ [" ... "] ++[last rr_]+           | otherwise = rr_+++toDenseListClip :: (Show a, Num a) => SpVector a -> Int -> String+toDenseListClip sv ncomax+  | dim sv > ncomax = unwords (map show h) ++  " ... " ++ show t+  | otherwise = show dr+     where dr = toDenseListSV sv+           h = take (ncomax - 2) dr+           t = last dr++printDenseSV :: (Show t, Num t) => SpVector t -> IO ()+printDenseSV sv = do+  newline+  printDenseSV' sv 5+  newline where+    printDenseSV' v nco = putStrLn rr_' where+      rr_ = toDenseListClip v nco :: String+      rr_' | dim sv > nco = unwords [take (nco - 2) rr_ , " ... " , [last rr_]]+           | otherwise = rr_++class PrintDense a where+  prd :: a -> IO ()++instance (Show a, Num a) => PrintDense (SpVector a) where+  prd = printDenseSV++instance (Show a, Num a) => PrintDense (SpMatrix a) where+  prd = printDenseSM++++++++++++++-- | rounding to 0 or 1 within some predefined numerical precision+++almostZero, almostOne :: Double -> Bool+almostZero x = abs x <= eps+almostOne x = x >= (1-eps) && x < (1+eps)++withDefault :: (t -> Bool) -> t -> t -> t+withDefault q d x | q x = d+                  | otherwise = x++roundZero, roundOne :: Double -> Double+roundZero = withDefault almostZero 0+roundOne = withDefault almostOne 1++with2Defaults :: (t -> Bool) -> (t -> Bool) -> t -> t -> t -> t+with2Defaults q1 q2 d1 d2 x | q1 x = d1+                            | q2 x = d2+                            | otherwise = x++roundZeroOne :: Double -> Double+roundZeroOne = with2Defaults almostZero almostOne 0 1++++++-- | transform state until a condition is met++modifyUntil :: MonadState s m => (s -> Bool) -> (s -> s) -> m s+modifyUntil q f = do+  x <- get+  let y = f x+  put y+  if q y then return y+         else modifyUntil q f     +  +loopUntilAcc :: Int -> ([t] -> Bool) -> (t -> t)  -> t -> t+loopUntilAcc nitermax q f x = go 0 [] x where+  go i ll xx | length ll < 2 = go (i + 1) (y : ll) y +             | otherwise = if q ll || i == nitermax+                           then xx+                           else go (i + 1) (take 2 $ y:ll) y+                where y = f xx++-- modify state and append, until max # of iterations is reached+modifyInspectN ::+  MonadState s m => Int -> ([s] -> Bool) -> (s -> s) -> m s+modifyInspectN nitermax q f +  | nitermax > 0 = go 0 []+  | otherwise = error "modifyInspectN : n must be > 0" where+      go i ll = do+        x <- get+        let y = f x+        if length ll < 2+          then do put y+                  go (i + 1) (y : ll)+          else if q ll || i == nitermax+               then do put y+                       return y+               else do put y+                       go (i + 1) (take 2 $ y : ll)+++meanl :: (Foldable t, Fractional a) => t a -> a+meanl xx = 1/fromIntegral (length xx) * sum xx++norm2l :: (Foldable t, Functor t, Floating a) => t a -> a+norm2l xx = sqrt $ sum (fmap (**2) xx)+++diffSqL xx = (x1 - x2)**2 where [x1, x2] = [head xx, xx!!1]+++++++++++untilConverged :: MonadState a m => (a -> SpVector Double) -> (a -> a) -> m a+untilConverged fproj = modifyInspectN 100 (normDiffConverged fproj)++-- convergence check (FIXME)+normDiffConverged :: (Foldable t, Functor t) =>+     (a -> SpVector Double) -> t a -> Bool+normDiffConverged fp xx = normSq (foldrMap fp (^-^) (zeroSV 0) xx) <= eps+              +++-- run `niter` iterations and append the state `x` to a list `xs`, stop when either the `xs` satisfies a predicate `q` or when the counter reaches 0++runAppendN :: ([t] -> Bool) -> (t -> t) -> Int -> t -> [t]+runAppendN qq ff niter x0 | niter<0 = error "runAppendN : niter must be > 0"+                          | otherwise = go qq ff niter x0 [] where+  go q f n z xs = +    let x = f z in+    if n <= 0 || q xs then xs+                      else go q f (n-1) x (x : xs)++-- ", NO convergence check +runAppendN' :: (t -> t) -> Int -> t -> [t]+runAppendN' ff niter x0 | niter<0 = error "runAppendN : niter must be > 0"+                        | otherwise = go ff niter x0 [] where+  go f n z xs = +    let x = f z in+    if n <= 0 then xs+              else go f (n-1) x (x : xs)++-- runN :: Int -> (a -> a) -> a -> a+-- runN n stepf x0 = runAppendN' stepf n x0+  ++++++++++-- | random matrices and vectors++-- dense++randMat :: PrimMonad m => Int -> m (SpMatrix Double)+randMat n = do+  g <- MWC.create+  aav <- replicateM (n^2) (MWC.normal 0 1 g)+  let ii_ = [0 .. n-1]+      (ix_,iy_) = unzip $ concatMap (zip ii_ . replicate n) ii_+  return $ fromListSM (n,n) $ zip3 ix_ iy_ aav+  +randVec :: PrimMonad m => Int -> m (SpVector Double)+randVec n = do+  g <- MWC.create+  bv <- replicateM n (MWC.normal 0 1 g)+  let ii_ = [0..n-1]+  return $ fromListSV n $ zip ii_ bv++++-- sparse++randSpMat :: Int -> Int -> IO (SpMatrix Double)+randSpMat n nsp | nsp > n^2 = error "randSpMat : nsp must be < n^2 "+                | otherwise = do+  g <- MWC.create+  aav <- replicateM nsp (MWC.normal 0 1 g)+  ii <- replicateM nsp (MWC.uniformR (0, n-1) g :: IO Int)+  jj <- replicateM nsp (MWC.uniformR (0, n-1) g :: IO Int)+  return $ fromListSM (n,n) $ zip3 ii jj aav+++randSpVec :: Int -> Int -> IO (SpVector Double)+randSpVec n nsp | nsp > n = error "randSpVec : nsp must be < n"+                | otherwise = do+  g <- MWC.create+  aav <- replicateM nsp (MWC.normal 0 1 g)+  ii <- replicateM nsp (MWC.uniformR (0, n-1) g :: IO Int)+  return $ fromListSV n $ zip ii aav+++++++-- | misc utils++++-- | integer-indexed ziplist+denseIxArray :: [b] -> [(Int, b)]+denseIxArray xs = zip [0..length xs-1] xs ++-- ", 2d arrays+denseIxArray2 :: Int -> [c] -> [(Int, Int, c)]+denseIxArray2 m xs = zip3 (concat $ replicate n ii_) jj_ xs where+  ii_ = [0 .. m-1]+  jj_ = concatMap (replicate m) [0 .. n-1]+  ln = length xs+  n = ln `div` m+++-- folds++foldrMap :: (Foldable t, Functor t) => (a -> b) -> (b -> c -> c) -> c -> t a -> c+foldrMap ff gg x0 = foldr gg x0 . fmap ff++foldlStrict :: (a -> b -> a) -> a -> [b] -> a+foldlStrict f = go+  where+    go z []     = z+    go z (x:xs) = let z' = f z x in z' `seq` go z' xs++-- indexed fold?+ifoldr :: Num i =>+     (a -> b -> b) -> b -> (i -> c -> d -> a) -> c -> [d] -> b  +ifoldr mjoin mneutral f  = go 0 where+  go i z (x:xs) = mjoin (f i z x) (go (i+1) z xs)+  go _ _ [] = mneutral+++-- bounds checking++type LB = Int+type UB = Int++inBounds :: LB -> UB -> Int -> Bool+inBounds ibl ibu i = i>= ibl && i<ibu++inBounds2 :: (LB, UB) -> (Int, Int) -> Bool+inBounds2 (ibl,ibu) (ix,iy) = inBounds ibl ibu ix && inBounds ibl ibu iy+++-- ", lower bound = 0+inBounds0 :: UB -> Int -> Bool+inBounds0 = inBounds 0++inBounds02 :: (UB, UB) -> (Int, Int) -> Bool+inBounds02 (bx,by) (i,j) = inBounds0 bx i && inBounds0 by j++++++++++--++tm0, tm1, tm2, tm3, tm4 :: SpMatrix Double+tm0 = fromListSM (2,2) [(0,0,pi), (1,0,sqrt 2), (0,1, exp 1), (1,1,sqrt 5)]++tv0, tv1 :: SpVector Double+tv0 = mkSpVectorD 2 [5, 6]+++tv1 = SV 2 $ IM.singleton 0 1++-- wikipedia test matrix for Givens rotation++tm1 = sparsifySM $ fromListDenseSM 3 [6,5,0,5,1,4,0,4,3]++tm1g1 = givens tm1 1 0+tm1a2 = tm1g1 ## tm1++tm1g2 = givens tm1a2 2 1+tm1a3 = tm1g2 ## tm1a2++tm1q = transposeSM (tm1g2 ## tm1g1)+++-- wp test matrix for QR decomposition via Givens rotation++tm2 = fromListDenseSM 3 [12, 6, -4, -51, 167, 24, 4, -68, -41]+++++tm3 = transposeSM $ fromListDenseSM 3 [1 .. 9]++tm3g1 = fromListDenseSM 3 [1, 0,0, 0,c,-s, 0, s, c]+  where c= 0.4961+        s = 0.8682+++--++tm4 = sparsifySM $ fromListDenseSM 4 [1,0,0,0,2,5,0,10,3,6,8,11,4,7,9,12]+++-- playground++-- | terminate after n iterations or when q becomes true, whichever comes first+untilC :: (a -> Bool) -> Int ->  (a -> a) -> a -> a+untilC p n f = go n+  where+    go m x | p x || m <= 0 = x+           | otherwise     = True `seq` go (m-1) (f x)++++++-- testing State+++-- data T0 = T0 {unT :: Int} deriving Eq+-- instance Show T0 where+--   show (T0 x) = show x++-- -- modifyT :: MonadState T0 m => (Int -> Int) -> m String+-- modifyT f = state (\(T0 i) -> (i, T0 (f i)))+  ++-- t00 = T0 0++-- testT n = execState $ replicateM n (modifyT (+1)) +++-- testT2 = execState $ when +  ++-- replicateSwitch p m f = loop m where+--       loop n | n <= 0 || p = pure (#)+--              | otherwise = f *> loop (n-1)++++-- testing Writer+               +-- asdfw n = runWriter $ do+--   tell $ "potato " ++ show n+--   tell "jam"+--   return (n+1)+++-- --+++-- testing State and Writer++++-- runMyApp runA k maxDepth =+--     let config = maxDepth+--         state =  0+--     in runStateT (runWriterT (runA k) config) state+++  

@@ -0,0 +1,149 @@+module Math.Linear.Sparse.IntMap where++import qualified Data.IntMap.Strict as IM+++++-- | ========= IntMap-of-IntMap (IM2) stuff+++-- insert an element+insertIM2 ::+  IM.Key -> IM.Key -> a -> IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a)+insertIM2 i j x imm = IM.insert i ro imm where+  ro = maybe (IM.singleton j x) (IM.insert j x) (IM.lookup i imm)+{-# inline insertIM2 #-}  ++-- lookup a key+lookupIM2 ::+  IM.Key -> IM.Key -> IM.IntMap (IM.IntMap a) -> Maybe a+lookupIM2 i j imm = IM.lookup i imm >>= IM.lookup j+{-# inline lookupIM2 #-}  ++-- populate an IM2 from a list of (row index, column index, value)  +fromListIM2 ::+  Foldable t =>+     t (IM.Key, IM.Key, a) -> IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a)+fromListIM2 iix sm = foldl ins sm iix where+  ins t (i,j,x) = insertIM2 i j x t+++-- | folding++-- indexed fold over an IM2+ifoldlIM2' :: (IM.Key -> IM.Key -> a -> b -> b) -> b -> IM.IntMap (IM.IntMap a) -> b+ifoldlIM2' f empty mm = IM.foldlWithKey' accRow empty mm where+  accRow acc i r = IM.foldlWithKey' (accElem i) acc r+  accElem i acc j x = f i j x acc+{-# inline ifoldlIM2' #-}++ifoldlIM2 ::+  (IM.Key -> IM.Key -> t -> IM.IntMap a -> IM.IntMap a) ->+  IM.IntMap (IM.IntMap t) ->  +  IM.IntMap a+ifoldlIM2 f m         = IM.foldlWithKey' accRow IM.empty m where+  accRow    acc i row = IM.foldlWithKey' (accElem i) acc row+  accElem i acc j x   = f i j x acc+{-# inline ifoldlIM2 #-}  ++foldlIM2 :: (a -> b -> b) -> b -> IM.IntMap (IM.IntMap a) -> b+foldlIM2 f empty mm = IM.foldl accRow empty mm where+  accRow acc r = IM.foldl accElem acc r+  accElem acc x = f x acc+{-# inline foldlIM2 #-}+++-- transposeIM2 : inner indices become outer ones and vice versa. No loss of information because both inner and outer IntMaps are nubbed.+transposeIM2 :: IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a)+transposeIM2 = ifoldlIM2 (flip insertIM2)+{-# inline transposeIM2 #-}++-- specialized folds++-- -- extract diagonal elements+-- extractDiagonalIM2 :: IM.IntMap (IM.IntMap a) -> [a]+-- extractDiagonalIM2 = ifoldlIM2' (\i j x xs -> if i==j then x : xs else xs) []+++++-- | filtering++-- map over outer IM and filter all inner IM's+ifilterIM2 ::+  (IM.Key -> IM.Key -> a -> Bool) ->+  IM.IntMap (IM.IntMap a) ->+  IM.IntMap (IM.IntMap a)+ifilterIM2 f  =+  IM.mapWithKey (\irow row -> IM.filterWithKey (f irow) row) +{-# inline ifilterIM2 #-}++-- specialized filtering function++-- keep only sub-diagonal elements+filterSubdiag :: IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a)+filterSubdiag = ifilterIM2 (\i j _ -> i>j)++countSubdiagonalNZ :: IM.IntMap (IM.IntMap a) -> Int+countSubdiagonalNZ im =+  IM.size $ IM.filter (not . IM.null) (filterSubdiag im)++-- list of (row, col) indices of (nonzero) subdiagonal elements+subdiagIndices :: IM.IntMap (IM.IntMap a) -> [(IM.Key, IM.Key)]+subdiagIndices im = concatMap rpairs $ IM.toList (IM.map IM.keys im') where+  im' = filterSubdiag im++rpairs :: (a, [b]) -> [(a, b)]+rpairs (i, jj@(_:_)) = zip (replicate (length jj) i) jj+rpairs (_, []) = []++-- -- list of (row, col) indices of elements that satisfy a criterion+-- indicesThatIM2 ::+--   (IM.Key -> IM.IntMap a -> Bool) -> IM.IntMap (IM.IntMap a) -> [(IM.Key, IM.Key)]+-- indicesThatIM2 f im = concatMap rpairs $ IM.toList (IM.map IM.keys im') where+--   im' = IM.filterWithKey f im++  +++-- | mapping++-- map over IM2++mapIM2 :: (a -> b) -> IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap b)+mapIM2 = IM.map . IM.map   -- imapIM2 (\_ _ x -> f x)+++-- indexed map over IM2+imapIM2 ::+  (IM.Key -> IM.Key -> a -> b) ->+  IM.IntMap (IM.IntMap a) ->+  IM.IntMap (IM.IntMap b)+imapIM2 f im = IM.mapWithKey ff im where+  ff j x = IM.mapWithKey (`f` j) x++++-- mapping keys++mapKeysIM2 ::+  (IM.Key -> IM.Key) -> (IM.Key -> IM.Key) -> IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a)+mapKeysIM2 fi fj im = IM.map adjCols adjRows where+  adjRows = IM.mapKeys fi im+  adjCols = IM.mapKeys fj +++++-- map over a single `column`++mapColumnIM2 :: (b -> b) -> IM.IntMap (IM.IntMap b) -> Int -> IM.IntMap (IM.IntMap b)+mapColumnIM2 f im jj = imapIM2 (\i j x -> if j == jj then f x else x) im++++++-- sparsification :+

@@ -0,0 +1,1 @@+{-# OPTIONS_GHC -F -pgmF hspec-discover #-}