diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..56eb9ea9bf7e2e450bb982763d4a88b1fce77fc7
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,46 @@
+*.py[co]
+
+# Packages
+*.egg
+*.egg-info
+dist
+build
+eggs
+parts
+bin
+var
+sdist
+develop-eggs
+.installed.cfg
+
+# Installer logs
+pip-log.txt
+
+# Mac OS X files
+.DS_Store
+.DS_Store?
+._*
+.Spotlight-V100
+.Trashes
+ehthumbs.db
+Thumbs.db
+
+# Sphinx stuff
+docs/generated/
+docs/build/
+
+# jupyter notebook
+.ipynb_checkpoints
+
+# coverage
+htmlcov
+.coverage
+
+# ide 
+.idea
+.spyderproject
+
+# Misc
+.cache
+.pytest_cache
+
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
new file mode 100644
index 0000000000000000000000000000000000000000..08e3844621f16b025513b84caa241df839a7f695
--- /dev/null
+++ b/.gitlab-ci.yml
@@ -0,0 +1,24 @@
+# run the test suite
+tests:
+    image: registry.gitlab.lis-lab.fr:5005/skmad-suite/madarrays/ubuntu:18.04
+    tags:
+        - docker
+    script:
+        - pip3 install --no-deps .
+        - pytest-3
+
+# generate the documentation
+pages:
+    image: registry.gitlab.lis-lab.fr:5005/skmad-suite/madarrays/ubuntu:18.04
+    tags:
+        - docker
+    only:
+        - master
+    script:
+        - pip3 install --no-deps .
+        - python3 setup.py build_sphinx
+        - cp -r build/sphinx/html public
+    artifacts:
+        paths:
+            - public
+
diff --git a/LICENSE.txt b/LICENSE.txt
new file mode 100755
index 0000000000000000000000000000000000000000..810fce6e9bf2aa10265b85614db5ac65941ecf81
--- /dev/null
+++ b/LICENSE.txt
@@ -0,0 +1,621 @@
+                    GNU GENERAL PUBLIC LICENSE
+                       Version 3, 29 June 2007
+
+ Copyright (C) 2007 Free Software Foundation, Inc. <https://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
diff --git a/MANIFEST.in b/MANIFEST.in
new file mode 100755
index 0000000000000000000000000000000000000000..eeb4dd5c0f6f9b91c8ceaa68179c33f9a8b31ec7
--- /dev/null
+++ b/MANIFEST.in
@@ -0,0 +1,8 @@
+include *.txt
+include *.rst
+include VERSION
+recursive-include doc *.rst *.py *.ipynb
+
+include pomad/tests/*.py
+
+prune doc/build
diff --git a/README.md b/README.md
deleted file mode 100644
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000
diff --git a/README.rst b/README.rst
new file mode 100755
index 0000000000000000000000000000000000000000..6e0a0308169022e93448395886450ed72767198b
--- /dev/null
+++ b/README.rst
@@ -0,0 +1,61 @@
+pomad
+=====
+
+pomad: PrObabilistic MAtrix Decompositions from Halko et al., 2011
+
+Python implementation of algorithms from paper *Finding Structure
+with Randomness: Probabilistic Algorithms for Constructing Approximate Matrix
+Decompositions*, by N. Halko, P. G. Martinsson and J. A. Tropp, SIAM review,
+53 (2), 2011.
+
+
+Install
+-------
+
+Install the current release with ``pip``::
+
+    pip install pomad
+
+For additional details, see doc/install.rst.
+
+Usage
+-----
+
+See the `documentation <http://templates.pages.lis-lab.fr/pomad/>`_.
+
+Bugs
+----
+
+Please report any bugs that you find through the `pomad GitLab project
+<https://gitlab.lis-lab.fr/valentin.emiya/pomad/issues>`_.
+
+You can also fork the repository and create a merge request.
+
+Source code
+-----------
+
+The source code of pomad is available via its `GitLab project
+<https://gitlab.lis-lab.fr/valentin.emiya/pomad>`_.
+
+You can clone the git repository of the project using the command::
+
+    git clone git@gitlab.lis-lab.fr:valentin.emiya/pomad.git
+
+Copyright © 2019-2020
+---------------------
+
+* `Laboratoire d'Informatique et Systèmes <http://www.lis-lab.fr/>`_
+* `Université d'Aix-Marseille <http://www.univ-amu.fr/>`_
+* `Centre National de la Recherche Scientifique <http://www.cnrs.fr/>`_
+* `Université de Toulon <http://www.univ-tln.fr/>`_
+
+Contributors
+------------
+
+* `Valentin Emiya <mailto:valentin.emiya@lis-lab.fr>`_
+
+License
+-------
+
+Released under the GNU General Public License version 3 or later
+(see `LICENSE.txt`).
diff --git a/VERSION b/VERSION
new file mode 100755
index 0000000000000000000000000000000000000000..6eec302096c96f3b65121bf64283ea6f63f11f22
--- /dev/null
+++ b/VERSION
@@ -0,0 +1 @@
+pomad:0.1
diff --git a/doc/Makefile b/doc/Makefile
new file mode 100755
index 0000000000000000000000000000000000000000..64227f569466b8601987972e617aabe461d39d71
--- /dev/null
+++ b/doc/Makefile
@@ -0,0 +1,121 @@
+# Makefile for Sphinx documentation
+#
+
+# You can set these variables from the command line.
+SPHINXOPTS    =
+SPHINXBUILD   = python -m sphinx 
+PAPER         =
+
+# Internal variables.
+PAPEROPT_a4     = -D latex_paper_size=a4
+PAPEROPT_letter = -D latex_paper_size=letter
+ALLSPHINXOPTS   = -d build/doctrees $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) .
+
+.PHONY: help clean html dirhtml pickle json htmlhelp qthelp latex changes linkcheck doctest epub
+
+help:
+	@echo "Please use \`make <target>' where <target> is one of"
+	@echo "  html      to make standalone HTML files"
+	@echo "  dirhtml   to make HTML files named index.html in directories"
+	@echo "  pickle    to make pickle files"
+	@echo "  epub       to make an epub"
+	@echo "  json      to make JSON files"
+	@echo "  htmlhelp  to make HTML files and a HTML help project"
+	@echo "  qthelp    to make HTML files and a qthelp project"
+	@echo "  latex     to make LaTeX files, you can set PAPER=a4 or PAPER=letter"
+	@echo "  changes   to make an overview of all changed/added/deprecated items"
+	@echo "  linkcheck to check all external links for integrity"
+	@echo "  doctest   to run all doctests embedded in the documentation (if enabled)"
+	@echo "  gitwash   to update the gitwash documentation"
+
+
+clean:
+	-rm -rf build/*
+	-rm -rf ghpages_build
+	-rm -rf auto_examples modules
+	-rm -rf reference/generated reference/algorithms/generated reference/classes/generated reference/readwrite/generated
+
+dist: html
+	test -d build/latex || make latex
+	make -C build/latex all-pdf
+	-rm -rf build/dist
+	(cd build/html; cp -r . ../../build/dist)
+	(cd build/dist && tar czf ../dist.tar.gz .)
+
+html:
+	$(SPHINXBUILD) -b html $(ALLSPHINXOPTS) build/html
+	@echo
+	@echo "Build finished. The HTML pages are in build/html."
+
+dirhtml:
+	$(SPHINXBUILD) -b dirhtml $(ALLSPHINXOPTS) build/dirhtml
+	@echo
+	@echo "Build finished. The HTML pages are in build/dirhtml."
+
+pickle:
+	$(SPHINXBUILD) -b pickle $(ALLSPHINXOPTS) build/pickle
+	@echo
+	@echo "Build finished; now you can process the pickle files."
+
+json:
+	$(SPHINXBUILD) -b json $(ALLSPHINXOPTS) build/json
+	@echo
+	@echo "Build finished; now you can process the JSON files."
+
+htmlhelp:
+	$(SPHINXBUILD) -b htmlhelp $(ALLSPHINXOPTS) build/htmlhelp
+	@echo
+	@echo "Build finished; now you can run HTML Help Workshop with the" \
+	      ".hhp project file in build/htmlhelp."
+
+qthelp:
+	$(SPHINXBUILD) -b qthelp $(ALLSPHINXOPTS) build/qthelp
+	@echo
+	@echo "Build finished; now you can run "qcollectiongenerator" with the" \
+	      ".qhcp project file in build/qthelp, like this:"
+	@echo "# qcollectiongenerator build/qthelp/test.qhcp"
+	@echo "To view the help file:"
+	@echo "# assistant -collectionFile build/qthelp/test.qhc"
+
+epub:
+	$(SPHINXBUILD) -b epub $(ALLSPHINXOPTS) build/epub
+	@echo
+	@echo "Build finished. The epub file is in $(BUILDDIR)/epub."
+
+
+latex:
+	$(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) build/latex
+	@echo
+	@echo "Build finished; the LaTeX files are in build/latex."
+	@echo "Run \`make all-pdf' or \`make all-ps' in that directory to" \
+	      "run these through (pdf)latex."
+
+changes:
+	$(SPHINXBUILD) -b changes $(ALLSPHINXOPTS) build/changes
+	@echo
+	@echo "The overview file is in build/changes."
+
+linkcheck:
+	$(SPHINXBUILD) -b linkcheck $(ALLSPHINXOPTS) build/linkcheck
+	@echo
+	@echo "Link check complete; look for any errors in the above output " \
+	      "or in build/linkcheck/output.txt."
+
+doctest:
+	$(SPHINXBUILD) -b doctest $(ALLSPHINXOPTS) build/doctest
+	@echo "Testing of doctests in the sources finished, look at the " \
+	      "results in build/doctest/output.txt."
+
+latexpdf: latex
+	@echo "Running LaTeX files through latexmk..."
+	$(MAKE) -C build/latex all-pdf
+	@echo "latexmk finished; the PDF files are in build/latex."
+
+docs: clean html latexpdf
+	cp build/latex/networkx_reference.pdf build/html/_downloads/.
+	
+gitwash-update:
+	python ../tools/gitwash_dumper.py developer networkx \
+	--project-url=http://networkx.github.io \
+	--project-ml-url=http://groups.google.com/group/networkx-discuss/ \
+	--gitwash-url git@github.com:matthew-brett/gitwash.git
diff --git a/doc/README.md b/doc/README.md
new file mode 100755
index 0000000000000000000000000000000000000000..5d2c9e95f0f138053f2d966c3d07e7de1029cab4
--- /dev/null
+++ b/doc/README.md
@@ -0,0 +1,28 @@
+If you only want to get the documentation, note that a pre-built
+version for the latest release is available
+[online](http://templates.pages.lis-lab.fr/pomad/).
+
+Sphinx is used to generate the API and reference documentation.
+
+## Instructions to build the documentation
+
+In addition to installing ``pomad`` and its dependencies, install the
+Python packages needed to build the documentation by entering
+
+```
+pip install -r ../requirements/doc.txt
+```
+in the ``doc/`` directory.
+
+To build the HTML documentation, run:
+```
+make html
+```
+in the ``doc/`` directory. This will generate a ``build/html`` subdirectory
+containing the built documentation.
+
+To build the PDF documentation, run:
+```
+make latexpdf
+```
+You will need to have Latex installed for this.
diff --git a/doc/_notebooks/pomad.ipynb b/doc/_notebooks/pomad.ipynb
new file mode 100755
index 0000000000000000000000000000000000000000..950b15439c8572bc2cc5a4358eeea204392c8d45
--- /dev/null
+++ b/doc/_notebooks/pomad.ipynb
@@ -0,0 +1,189 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Tutorial for `pomad`:  PrObabilistic MAtrix Decompositions"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "%load_ext autoreload\n",
+    "%autoreload 2"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "%matplotlib inline"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "import scipy as sp\n",
+    "from scipy.sparse.linalg import aslinearoperator"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "from pomad.utils import build_test_matrix\n",
+    "from pomad.range_approximation import randomized_range_finder, adaptive_randomized_range_finder"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "# Test behavior of build_test_matrix\n",
+    "m = 29\n",
+    "n = 17\n",
+    "seed = 0\n",
+    "p_max = 5\n",
+    "rand_state = np.random.RandomState(seed)\n",
+    "A = rand_state.randn(m, n)\n",
+    "A = A / np.linalg.norm(A, ord=2)\n",
+    "M = A.copy()\n",
+    "for p in range(p_max):\n",
+    "    U, s, Vh = np.linalg.svd(M)\n",
+    "    plt.semilogy(s, label=r'$(A A^T)^{} \\times A$'.format(p))\n",
+    "    M = A @ A.T @ M\n",
+    "plt.title(r'Singular values of $(A A^T)^p \\times A$')\n",
+    "plt.legend()\n",
+    "plt.grid()\n",
+    "\n",
+    "plt.figure()\n",
+    "for p in range(p_max):\n",
+    "    M = build_test_matrix(m, n, p=p, rand_state=seed)\n",
+    "    U, s, Vh = np.linalg.svd(M)\n",
+    "    plt.semilogy(s, label='p={}'.format(p))\n",
+    "plt.title('Similar result using build_test_matrix(m, n, p)')\n",
+    "plt.legend()\n",
+    "plt.grid()\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "# Approximation quality of randomized_range_finder(A, n_l=l)\n",
+    "m = 29\n",
+    "n = 23\n",
+    "p = 3\n",
+    "A = build_test_matrix(m, n, p=10, rand_state=0)\n",
+    "_, s, _ = np.linalg.svd(A)\n",
+    "for k in range(1, 18):\n",
+    "    Q = randomized_range_finder(A, n_l=k+p)\n",
+    "    print(\"k =\", k, \"p =\", p, \n",
+    "          \"||A - Q @ Q.T @ A|| < s[k]?\", np.linalg.norm(A - Q @ Q.T @ A, ord=2) < s[k], \n",
+    "          \"({} < {})\".format(np.linalg.norm(A - Q @ Q.T @ A, ord=2), s[k]))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "m = 290\n",
+    "n = 170\n",
+    "A = build_test_matrix(m, n, p=10, rand_state=0)\n",
+    "_, s, _ = np.linalg.svd(A)\n",
+    "plt.figure(figsize=(15,5))\n",
+    "plt.semilogy(s)\n",
+    "plt.grid()\n",
+    "for l in range(10):\n",
+    "    tolerance = 10**-l\n",
+    "    Q = adaptive_randomized_range_finder(A, tolerance=tolerance, r=5)\n",
+    "    print(tolerance, np.linalg.norm(A-Q@Q.T@A, ord=2), Q.shape[1])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "m = n = 345\n",
+    "A = build_test_matrix(m, n, p=10, rand_state=0)\n",
+    "Aop = aslinearoperator(A)\n",
+    "tolerance = 10**-5\n",
+    "%timeit -r2 -p2 Q = adaptive_randomized_range_finder(Aop, tolerance=tolerance, r=5)\n",
+    "A_rank = Q.shape[1]\n",
+    "print(A_rank)\n",
+    "%timeit -r2 -p2 QQ = randomized_range_finder(Aop, n_l=A_rank)\n",
+    "QQ = randomized_range_finder(Aop, n_l=A_rank)\n",
+    "%timeit -r2 -p2 np.linalg.svd(np.conjugate(QQ).T @ A)\n",
+    "%timeit -r2 -p2 ss=sp.sparse.linalg.svds(Aop, k=A_rank)\n",
+    "%timeit -r2 -p2 sss=np.linalg.svd(A)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.6.2"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/doc/conf.py b/doc/conf.py
new file mode 100755
index 0000000000000000000000000000000000000000..addf02bcb697a696b567ecc3206e9c24b78a5a7c
--- /dev/null
+++ b/doc/conf.py
@@ -0,0 +1,286 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+import sys
+import os
+
+from datetime import date
+
+import pomad
+
+# If extensions (or modules to document with autodoc) are in another directory,
+# add these directories to sys.path here. If the directory is relative to the
+# documentation root, use os.path.abspath to make it absolute, like shown here.
+sys.path.insert(0, os.path.abspath('../pomad/'))
+
+# -- General configuration ------------------------------------------------
+
+# Add any Sphinx extension module names here, as strings. They can be
+# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom ones.
+
+extensions = ['sphinx.ext.autodoc',
+              'sphinx.ext.autosummary',
+              'sphinx.ext.doctest',
+              'sphinx.ext.intersphinx',
+              'sphinx.ext.todo',
+              'sphinx.ext.coverage',
+              'sphinx.ext.mathjax',
+              'sphinx.ext.viewcode',
+              'numpydoc',
+              'nbsphinx',
+              'IPython.sphinxext.ipython_console_highlighting',
+              ]
+
+
+# Add any paths that contain templates here, relative to this directory.
+templates_path = ['_templates']
+
+# The suffix(es) of source filenames.
+# You can specify multiple suffix as a list of string:
+# source_suffix = ['.rst', '.md']
+source_suffix = '.rst'
+
+# The encoding of source files.
+source_encoding = 'utf-8-sig'
+
+# The master toctree document.
+master_doc = 'index'
+
+# General information about the project.
+project = 'pomad'
+author = 'V. Emiya'
+copyright = '2019-{}, {}'.format(date.today().year, author)
+
+# The version info for the project you're documenting, acts as replacement for
+# |version| and |release|, also used in various other places throughout the
+# built documents.
+#
+# The short X.Y version.
+version = pomad.__version__
+# The full version, including alpha/beta/rc tags.
+release = pomad.__version__.replace('_', '')
+
+# The language for content autogenerated by Sphinx. Refer to documentation
+# for a list of supported languages.
+#
+# This is also used if you do content translation via gettext catalogs.
+# Usually you set "language" from the command line for these cases.
+language = None
+
+# There are two options for replacing |today|: either, you set today to some
+# non-false value, then it is used:
+# Else, today_fmt is used as the format for a strftime call.
+today_fmt = '%B %d, %Y'
+
+# List of patterns, relative to source directory, that match files and
+# directories to ignore when looking for source files.
+exclude_patterns = ['_build', '**/test_*.rst', '**.ipynb_checkpoints']
+# The reST default role (used for this markup: `text`) to use for all
+# documents.
+default_role = None
+
+# If true, '()' will be appended to :func: etc. cross-reference text.
+add_function_parentheses = True
+
+# If true, the current module name will be prepended to all description
+# unit titles (such as .. function::).
+add_module_names = True
+
+# If true, sectionauthor and moduleauthor directives will be shown in the
+# output. They are ignored by default.
+show_authors = False
+
+# The name of the Pygments (syntax highlighting) style to use.
+pygments_style = 'sphinx'
+
+# A list of ignored prefixes for module index sorting.
+modindex_common_prefix = []
+
+# If true, keep warnings as "system message" paragraphs in the built documents.
+keep_warnings = False
+
+# If true, `todo` and `todoList` produce output, else they produce nothing.
+todo_include_todos = True
+
+# -- Options for HTML output ----------------------------------------------
+
+# The theme to use for HTML and HTML Help pages.  See the documentation for
+# a list of builtin themes.
+html_theme = 'bizstyle'
+
+# Theme options are theme-specific and customize the look and feel of a theme
+# further.  For a list of options available for each theme, see the
+# documentation.
+html_theme_options = {}
+
+# Add any paths that contain custom themes here, relative to this directory.
+html_theme_path = []
+
+# The name for this set of Sphinx documents.  If None, it defaults to
+# "<project> v<release> documentation".
+# html_title = None
+
+# A shorter title for the navigation bar.  Default is the same as html_title.
+# html_short_title = None
+
+# The name of an image file (relative to this directory) to place at the top
+# of the sidebar.
+# html_logo = None
+
+# The name of an image file (within the static path) to use as favicon of the
+# docs.  This file should be a Windows icon file (.ico) being 16x16 or 32x32
+# pixels large.
+html_favicon = None
+
+# Add any paths that contain custom static files (such as style sheets) here,
+# relative to this directory. They are copied after the builtin static files,
+# so a file named "default.css" will overwrite the builtin "default.css".
+html_static_path = []
+
+# Add any extra paths that contain custom files (such as robots.txt or
+# .htaccess) here, relative to this directory. These files are copied
+# directly to the root of the documentation.
+html_extra_path = []
+
+# If not '', a 'Last updated on:' timestamp is inserted at every page bottom,
+# using the given strftime format.
+html_last_updated_fmt = '%b %d, %Y'
+
+# Custom sidebar templates, maps document names to template names.
+html_sidebars = {}
+
+# Additional templates that should be rendered to pages, maps page names to
+# template names.
+html_additional_pages = {}
+
+# If false, no module index is generated.
+html_domain_indices = True
+
+# If false, no index is generated.
+html_use_index = True
+
+# If true, the index is split into individual pages for each letter.
+html_split_index = False
+
+# If true, links to the reST sources are added to the pages.
+html_show_sourcelink = True
+
+# If true, "Created using Sphinx" is shown in the HTML footer. Default is True.
+html_show_sphinx = True
+
+# If true, "(C) Copyright ..." is shown in the HTML footer. Default is True.
+html_show_copyright = True
+
+# If true, an OpenSearch description file will be output, and all pages will
+# contain a <link> tag referring to it.  The value of this option must be the
+# base URL from which the finished HTML is served.
+html_use_opensearch = ''
+
+# This is the file name suffix for HTML files (e.g. ".xhtml").
+html_file_suffix = None
+
+# Language to be used for generating the HTML full-text search index.
+# Sphinx supports the following languages:
+#   'da', 'de', 'en', 'es', 'fi', 'fr', 'h', 'it', 'ja'
+#   'nl', 'no', 'pt', 'ro', 'r', 'sv', 'tr'
+html_search_language = 'en'
+
+# A dictionary with options for the search language support, empty by default.
+# Now only 'ja' uses this config value
+html_search_options = {'type': 'default'}
+
+# The name of a javascript file (relative to the configuration directory) that
+# implements a search results scorer. If empty, the default will be used.
+html_search_scorer = ''
+
+# Output file base name for HTML help builder.
+htmlhelp_basename = 'pomaddoc'
+
+# -- Options for LaTeX output ---------------------------------------------
+
+latex_elements = {
+    # The paper size ('letterpaper' or 'a4paper').
+    'papersize': 'a4paper',
+
+    # The font size ('10pt', '11pt' or '12pt').
+    'pointsize': '10pt',
+
+    # Additional stuff for the LaTeX preamble.
+    'preamble': '',
+
+    # Latex figure (float) alignment
+    'figure_align': 'htbp',
+}
+
+# Grouping the document tree into LaTeX files. List of tuples
+# (source start file, target name, title,
+#  author, documentclass [howto, manual, or own class]).
+latex_documents = [
+    (master_doc, 'pomad.tex', 'pomad Documentation', author, 'manual'),
+]
+
+# The name of an image file (relative to this directory) to place at the top of
+# the title page.
+latex_logo = None
+
+# For "manual" documents, if this is true, then toplevel headings are parts,
+# not chapters.
+latex_use_parts = False
+
+# If true, show page references after internal links.
+latex_show_pagerefs = False
+
+# If true, show URL addresses after external links.
+# latex_show_urls = False
+
+# Documents to append as an appendix to all manuals.
+latex_appendices = []
+
+# If false, no module index is generated.
+latex_domain_indices = True
+
+# -- Options for manual page output ---------------------------------------
+
+# One entry per manual page. List of tuples
+# (source start file, name, description, authors, manual section).
+man_pages = [(master_doc, 'pomad', 'pomad Documentation', [author], 1)]
+
+# If true, show URL addresses after external links.
+man_show_urls = False
+
+# -- Options for Texinfo output -------------------------------------------
+
+# Grouping the document tree into Texinfo files. List of tuples
+# (source start file, target name, title, author,
+#  dir menu entry, description, category)
+texinfo_documents = [
+    (master_doc, 'pomad.tex', 'pomad Documentation', author, 'manual'),
+]
+
+# Documents to append as an appendix to all manuals.
+texinfo_appendices = []
+
+# If false, no module index is generated.
+texinfo_domain_indices = True
+
+# How to display URL addresses: 'footnote', 'no', or 'inline'.
+texinfo_show_urls = 'footnote'
+
+# If true, do not generate a @detailmenu in the "Top" node's menu.
+texinfo_no_detailmenu = False
+
+# Example configuration for intersphinx: refer to the Python standard library.
+intersphinx_mapping = {
+    'numpy': ('https://docs.scipy.org/doc/numpy/', None),
+}
+
+# Allow errors in notebook
+nbsphinx_allow_errors = True
+
+# Do not show class members
+numpydoc_show_class_members = False
+
+# Include todos
+todo_include_todos = True
+
+# Order members by member type
+autodoc_member_order = 'groupwise'
diff --git a/doc/copyright.rst b/doc/copyright.rst
new file mode 100755
index 0000000000000000000000000000000000000000..180e86a44689775ac74ed164e0bcb663615ca5b7
--- /dev/null
+++ b/doc/copyright.rst
@@ -0,0 +1,49 @@
+Credits
+#######
+
+Copyright(c) 2019-2020
+----------------------
+
+* Laboratoire d'Informatique et Systèmes <http://www.lis-lab.fr/>
+* Université d'Aix-Marseille <http://www.univ-amu.fr/>
+* Centre National de la Recherche Scientifique <http://www.cnrs.fr/>
+* Université de Toulon <http://www.univ-tln.fr/>
+
+Contributors
+------------
+
+* Valentin Emiya <firstname.lastname_AT_lis-lab.fr>
+
+This package has been created thanks to the joint work with Florent Jaillet
+and Ronan Hamon on other packages.
+
+Description
+-----------
+
+`pomad` is a Python implementation of algorithms from
+paper *Finding Structure with Randomness: Probabilistic Algorithms for
+Constructing Approximate Matrix Decompositions*, by N. Halko, P. G.
+Martinsson and J. A. Tropp, SIAM review, 53 (2), 2011, https://arxiv.org/abs/0909.4061.
+
+
+Version
+-------
+
+* pomad version = 0.1
+
+Licence
+-------
+
+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/>.
+
diff --git a/doc/credits.rst b/doc/credits.rst
new file mode 100755
index 0000000000000000000000000000000000000000..80f497d5223e00a394cd9ba94530666d9770b868
--- /dev/null
+++ b/doc/credits.rst
@@ -0,0 +1,35 @@
+Credits
+=======
+
+Copyright(c) 2020
+-----------------
+
+* Laboratoire d'Informatique et Systèmes <http://www.lis-lab.fr/>
+* Université d'Aix-Marseille <http://www.univ-amu.fr/>
+* Centre National de la Recherche Scientifique <http://www.cnrs.fr/>
+* Université de Toulon <http://www.univ-tln.fr/>
+
+Contributors
+------------
+
+* Valentin Emiya <firstname.lastname_AT_lis-lab.fr>
+
+This package has been created thanks to the joint work with Florent Jaillet
+and Ronan Hamon on other packages.
+
+Licence
+-------
+This file is part of pomad.
+
+pomad 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/>.
diff --git a/doc/index.rst b/doc/index.rst
new file mode 100755
index 0000000000000000000000000000000000000000..9e42c85078a998c4993d055773df6544f5f09446
--- /dev/null
+++ b/doc/index.rst
@@ -0,0 +1,43 @@
+##########################
+:mod:`pomad` documentation
+##########################
+
+Overview
+========
+:mod:`pomad`: stands for PrObabilistic MAtrix Decompositions.
+
+Package :mod:`pomad` offers a Python implementation of algorithms from
+paper *Finding Structure with Randomness: Probabilistic Algorithms for
+Constructing Approximate Matrix Decompositions*, by N. Halko, P. G.
+Martinsson and J. A. Tropp, SIAM review, 53 (2), 2011, https://arxiv.org/abs/0909.4061.
+
+:mod:`pomad`: is mainly composed of two sub-modules:
+
+* :mod:`pomad.range_approximation` contains algorithms for Stage A Randomized Schemes for Approximating the Range
+* :mod:`pomad.factorization_construction` contains algorithms for Stage B Construction of Standard Factorizations
+
+Not all algorithms have been implemented yet.
+
+Documentation
+=============
+
+.. only:: html
+
+    :Release: |version|
+    :Date: |today|
+
+.. toctree::
+    :maxdepth: 1
+
+    installation
+    references
+    tutorials
+    copyright
+
+
+Indices and tables
+==================
+
+* :ref:`genindex`
+* :ref:`modindex`
+* :ref:`search`
diff --git a/doc/installation.rst b/doc/installation.rst
new file mode 100755
index 0000000000000000000000000000000000000000..24d8496b2ff78418e5289b720ca992c3c7cb5656
--- /dev/null
+++ b/doc/installation.rst
@@ -0,0 +1,75 @@
+Installation
+############
+
+``pomad`` requires the following packages:
+
+* `python >= 3.5 <https://wiki.python.org/moin/BeginnersGuide/Download>`_
+* `numpy >= 1.13 <http://www.numpy.org>`_
+
+Make sure your Python environment is properly configured. It is recommended to
+install ``pomad`` in a virtual environment.
+
+Release version
+---------------
+
+First, make sure you have the latest version of pip (the Python package
+manager) installed. If you do not, refer to the `Pip documentation
+<https://pip.pypa.io/en/stable/installing/>`_ and install ``pip`` first.
+
+Install the current release with ``pip``::
+
+    pip install pomad
+
+To upgrade to a newer release use the ``--upgrade`` flag::
+
+    pip install --upgrade pomad
+
+If you do not have permission to install software systemwide, you can install
+into your user directory using the ``--user`` flag::
+
+    pip install --user pomad
+
+Alternatively, you can manually download ``pomad`` from its `GitLab project
+<https://gitlab.lis-lab.fr/valentin.emiya/pomad>`_  or `PyPI
+<https://pypi.python.org/pypi/pomad>`_.  To install one of these versions,
+unpack it and run the following from the top-level source directory using the
+Terminal::
+
+    pip install .
+
+Development version
+-------------------
+
+If you have `Git <https://git-scm.com/>`_ installed on your system, it is also
+possible to install the development version of ``pomad``.
+
+Before installing the development version, you may need to uninstall the
+standard version of ``pomad`` using ``pip``::
+
+    pip uninstall pomad
+
+Clone the Git repository::
+
+    git clone git@gitlab.lis-lab.fr:valentin.emiya/pomad.git
+    cd pomad
+
+You may also need to install required packages::
+
+    pip install -r requirements/defaults.txt
+
+Then execute ``pip`` with flag ``-e`` to follow the development branch::
+
+    pip install -e .
+
+To update ``pomad`` at any time, in the same directory do::
+
+    git pull
+
+To run unitary tests, first install required packages::
+
+    pip install -r requirements/dev.txt
+
+and execute ``pytest``::
+
+    pytest
+
diff --git a/doc/references.rst b/doc/references.rst
new file mode 100755
index 0000000000000000000000000000000000000000..e3d591d684ed064ef03d4c5ca521ae95e6d74fb0
--- /dev/null
+++ b/doc/references.rst
@@ -0,0 +1,29 @@
+References
+==========
+
+    :Release: |release|
+    :Date: |today|
+
+pomad\.factorization_construction module
+----------------------------------------
+
+.. automodule:: pomad.factorization_construction
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+pomad\.range_approximation module
+---------------------------------
+
+.. automodule:: pomad.range_approximation
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+pomad\.utils module
+-------------------
+
+.. automodule:: pomad.utils
+    :members:
+    :undoc-members:
+    :show-inheritance:
diff --git a/doc/tutorials.rst b/doc/tutorials.rst
new file mode 100755
index 0000000000000000000000000000000000000000..273e21d8106adde60380afb8aa378f5d475cefb5
--- /dev/null
+++ b/doc/tutorials.rst
@@ -0,0 +1,7 @@
+Tutorials
+#########
+
+.. toctree::
+    :maxdepth: 1
+
+    _notebooks/pomad.ipynb
diff --git a/pomad/__init__.py b/pomad/__init__.py
new file mode 100755
index 0000000000000000000000000000000000000000..985033fc17438fadcb5362afbc4073e37ace8b44
--- /dev/null
+++ b/pomad/__init__.py
@@ -0,0 +1,14 @@
+# -*- coding: utf-8 -*-
+"""
+.. moduleauthor:: Valentin Emiya
+"""
+from .factorization_construction import direct_svd, evd_nystrom
+from .range_approximation import \
+    randomized_range_finder, adaptive_randomized_range_finder
+
+# TODO replace RandomState by BitGenerator
+
+__all__ = ['direct_svd', 'evd_nystrom',
+           'randomized_range_finder', 'adaptive_randomized_range_finder']
+
+__version__ = "0.1"
diff --git a/pomad/_dev_range_approximation.py b/pomad/_dev_range_approximation.py
new file mode 100755
index 0000000000000000000000000000000000000000..ee46029832d589f6e622241cf351c0d5fc60232e
--- /dev/null
+++ b/pomad/_dev_range_approximation.py
@@ -0,0 +1,348 @@
+# -*- coding: utf-8 -*-
+# ######### COPYRIGHT #########
+# Credits
+# #######
+#
+# Copyright(c) 2019-2020
+# ----------------------
+#
+# * Laboratoire d'Informatique et Systèmes <http://www.lis-lab.fr/>
+# * Université d'Aix-Marseille <http://www.univ-amu.fr/>
+# * Centre National de la Recherche Scientifique <http://www.cnrs.fr/>
+# * Université de Toulon <http://www.univ-tln.fr/>
+#
+# Contributors
+# ------------
+#
+# * Valentin Emiya <firstname.lastname_AT_lis-lab.fr>
+#
+# This package has been created thanks to the joint work with Florent Jaillet
+# and Ronan Hamon on other packages.
+#
+# Description
+# -----------
+#
+# `pomad` is a Python implementation of algorithms from
+# paper *Finding Structure with Randomness: Probabilistic Algorithms for
+# Constructing Approximate Matrix Decompositions*, by N. Halko, P. G.
+# Martinsson and J. A. Tropp, SIAM review, 53 (2), 2011, https://arxiv.org/abs/0909.4061.
+#
+#
+# Version
+# -------
+#
+# * pomad version = 0.1
+#
+# Licence
+# -------
+#
+# 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/>.
+#
+# ######### COPYRIGHT #########
+"""
+
+.. moduleauthor:: Valentin Emiya
+"""
+import numpy as np
+from scipy.sparse.linalg import aslinearoperator
+
+
+def adaptive_randomized_range_finder_naive(a, tolerance, r=5, rand_state=None):
+    """ Naive implementation for Algo 4.2 """
+    if rand_state is None:
+        rand_state = np.random.RandomState(None)
+    if np.issubdtype(type(rand_state), np.dtype(int).type):
+        rand_state = np.random.RandomState(rand_state)
+    m, n = a.shape
+    W = rand_state.randn(n, r)
+    Y = [a @ W[:, i] for i in range(r)]
+    j = 0
+    Q = np.empty((m, 0), dtype=a.dtype)
+    y_norms = [np.linalg.norm(Y[i]) for i in range(r)]
+
+    while np.max(y_norms[j:j + r]) > tolerance / (10 * np.sqrt(2 / np.pi)):
+        j += 1
+        Y[j-1] -= Q @ (np.conj(Q).T @ Y[j-1])  # Fixed bug with index j
+        q = Y[j-1] / np.linalg.norm(Y[j-1])  # Fixed bug with index j
+        Q = np.concatenate((Q, q[:, None]), axis=1)
+        w = rand_state.randn(n)
+        Aw = a @ w
+        y = Aw - Q @ (Q.T.conj() @ Aw)
+        Y.append(y)
+        y_norms.append(np.linalg.norm(y))
+        for i in range(j, j + r - 1):  # Fixed bug with final index
+            Y[i] -= q * np.dot(q.conj(), Y[i])
+    return Q
+
+
+def adaptive_randomized_range_finder_nolist(a, tolerance, r=5, rand_state=None):
+    """ Intermediate implementation for Algo 4.2:
+    adaptive_randomized_range_finder_naive enhanced using nd-array instead
+    of lists  """
+    if rand_state is None:
+        rand_state = np.random.RandomState(None)
+    if np.issubdtype(type(rand_state), np.dtype(int).type):
+        rand_state = np.random.RandomState(rand_state)
+    m, n = a.shape
+    W = rand_state.randn(n, r)
+    Y = a @ W
+    y_norms = np.linalg.norm(Y, axis=0)
+    j = 0
+    Q = np.empty((m, 0), dtype=a.dtype)
+
+    while np.max(y_norms[j:j + r]) > tolerance / (10 * np.sqrt(2 / np.pi)):
+        Y[:, j] -= Q @ (np.conj(Q).T @ Y[:, j])
+        q = Y[:, j] / np.linalg.norm(Y[:, j])
+        j += 1
+        Q = np.concatenate((Q, q[:, None]), axis=1)
+        w = rand_state.randn(n)
+        Aw = a @ w
+        y = Aw - Q @ (Q.T.conj() @ Aw)
+        Y = np.concatenate((Y, y[:, None]), axis=1)
+        y_norms = np.concatenate((y_norms, [np.linalg.norm(y)]))
+        for i in range(j, j + r - 1):
+            Y[:, i] -= q * np.dot(q.conj(), Y[:, i])
+    return Q
+
+
+def adaptive_randomized_range_finder_circy(a, tolerance, r=5, rand_state=None):
+    """ Intermediate implementation for Algo 4.2:
+    adaptive_randomized_range_finder_nolist enhanced using circular array
+    for Y  """
+    if rand_state is None:
+        rand_state = np.random.RandomState(None)
+    if np.issubdtype(type(rand_state), np.dtype(int).type):
+        rand_state = np.random.RandomState(rand_state)
+    m, n = a.shape
+    W = rand_state.randn(n, r)
+    Y = a @ W
+    y_norms = np.linalg.norm(Y, axis=0)
+    j = 0
+    j_incr = 0
+    Q = np.empty((m, 0), dtype=a.dtype)
+
+    while np.max(y_norms[j_incr:j_incr + r]) > tolerance / (10 * np.sqrt(2 / np.pi)):
+        j_incr += 1
+        Y[:, j] -= Q @ (np.conj(Q).T @ Y[:, j])
+        q = Y[:, j] / np.linalg.norm(Y[:, j])
+        j = (j + 1) % r
+        Q = np.concatenate((Q, q[:, None]), axis=1)
+        w = rand_state.randn(n)
+        A_w = a @ w
+        y = A_w - Q @ (Q.T.conj() @ A_w)
+        Y[:, j - 1] = y
+        y_norms = np.concatenate((y_norms, [np.linalg.norm(y)]))
+        for i in np.arange(j, j + r - 1) % r:
+            Y[:, i] -= q * np.dot(q.conj(), Y[:, i])
+    return Q
+
+
+def adaptive_randomized_range_finder_circynorms(a, tolerance, r=5,
+                                                rand_state=None):
+    """ Intermediate implementation for Algo 4.2:
+    adaptive_randomized_range_finder_circy enhanced using circular array
+    for y_norms  """
+    if rand_state is None:
+        rand_state = np.random.RandomState(None)
+    if np.issubdtype(type(rand_state), np.dtype(int).type):
+        rand_state = np.random.RandomState(rand_state)
+    m, n = a.shape
+    W = rand_state.randn(n, r)
+    Y = a @ W
+    y_norms = np.linalg.norm(Y, axis=0)
+    j = 0
+    Q = np.empty((m, 0), dtype=a.dtype)
+
+    while np.max(y_norms) > tolerance / (10 * np.sqrt(2 / np.pi)):
+        Y[:, j] -= Q @ (np.conj(Q).T @ Y[:, j])
+        q = Y[:, j] / np.linalg.norm(Y[:, j])
+        j = (j + 1) % r
+        Q = np.concatenate((Q, q[:, None]), axis=1)
+        w = rand_state.randn(n)
+        A_w = a @ w
+        y = A_w - Q @ (Q.T.conj() @ A_w)
+        Y[:, j - 1] = y
+        y_norms[j - 1] = np.linalg.norm(y)
+        for i in np.arange(j, j + r - 1) % r:
+            Y[:, i] -= q * np.dot(q.conj(), Y[:, i])
+    return Q
+
+
+def adaptive_randomized_range_finder_rmloop(a, tolerance, r=5,
+                                            rand_state=None):
+    """ Intermediate implementation for Algo 4.2:
+    adaptive_randomized_range_finder_circynorms enhanced removing internal
+    loop  """
+    if rand_state is None:
+        rand_state = np.random.RandomState(None)
+    if np.issubdtype(type(rand_state), np.dtype(int).type):
+        rand_state = np.random.RandomState(rand_state)
+    m, n = a.shape
+    W = rand_state.randn(n, r)
+    Y = a @ W
+    y_norms = np.linalg.norm(Y, axis=0)
+    j = 0
+    Q = np.empty((m, 0), dtype=a.dtype)
+
+    while np.max(y_norms) > tolerance / (10 * np.sqrt(2 / np.pi)):
+        Y[:, j] -= Q @ (Q.T.conj() @ Y[:, j])
+        q = Y[:, j] / np.linalg.norm(Y[:, j])
+        j = (j + 1) % r
+        Q = np.concatenate((Q, q[:, None]), axis=1)
+        w = rand_state.randn(n)
+        A_w = a @ w
+        Y -= np.outer(q, q.conj() @ Y)
+        y = A_w - Q @ (Q.T.conj() @ A_w)
+        Y[:, j - 1] = y
+        y_norms[j - 1] = np.linalg.norm(y)
+    return Q
+
+
+def adaptive_randomized_range_finder_rmconcQ(a, tolerance, r=5,
+                                             rand_state=None):
+    """ Intermediate implementation for Algo 4.2:
+    adaptive_randomized_range_finder_rmloop enhanced removing concatenation
+    for Q, with maximal size for Q (needs to be reduced)  """
+    if rand_state is None:
+        rand_state = np.random.RandomState(None)
+    if np.issubdtype(type(rand_state), np.dtype(int).type):
+        rand_state = np.random.RandomState(rand_state)
+    m, n = a.shape
+    W = rand_state.randn(n, r)
+    Y = a @ W
+    y_norms = np.linalg.norm(Y, axis=0)
+    j = 0
+    Q_width = 0
+    Q = np.empty((m, min(m, n)), dtype=a.dtype)
+
+    while np.max(y_norms) > tolerance / (10 * np.sqrt(2 / np.pi)):
+        if Q_width >= n:
+            # TODO to be tested
+            raise RuntimeError('Required tolerance not reached with maximum '
+                               'number of columns ({}). Singular '
+                               'values may decrease too slowly.'.format(n))
+        Y[:, j] -= Q[:, :Q_width] @ (Q[:, :Q_width].T.conj() @ Y[:, j])
+        q = Y[:, j] / np.linalg.norm(Y[:, j])
+        j = (j + 1) % r
+        Q_width += 1
+        Q[:, Q_width-1] = q
+        w = rand_state.randn(n)
+        A_w = a @ w
+        # Y -= np.outer(q, np.conj(np.conj(q) @ Y))
+        Y -= np.outer(q, q.conj() @ Y)
+        y = A_w - Q[:, :Q_width] @ (Q[:, :Q_width].T.conj() @ A_w)
+        Y[:, j - 1] = y
+        y_norms[j - 1] = np.linalg.norm(y)
+    return Q[:, :Q_width]
+
+
+def adaptive_randomized_range_finder_mem_alloc(a, tolerance, r=5,
+                                               rand_state=None, n_cols_Q=16):
+    """ Intermediate implementation for Algo 4.2:
+    adaptive_randomized_range_finder_rmconcQ enhanced using memory
+    allocation (subsequently optimized using resize in final version
+    adaptive_randomized_range_finder)
+    """
+    if rand_state is None:
+        rand_state = np.random.RandomState(None)
+    if np.issubdtype(type(rand_state), np.dtype(int).type):
+        rand_state = np.random.RandomState(rand_state)
+    if isinstance(a, np.ndarray):
+        a = aslinearoperator(a)
+    m, n = a.shape
+    W = rand_state.randn(n, r)
+    Y = a @ W
+    y_norms = np.linalg.norm(Y, axis=0)
+    j = 0
+    Q_width = 0
+    Q = np.empty((m, n_cols_Q), order='F', dtype=a.dtype)
+
+    while np.max(y_norms) > tolerance / (10 * np.sqrt(2 / np.pi)):
+        if Q_width >= n:
+            # TODO to be tested
+            raise RuntimeError('Required tolerance not reached with maximum '
+                               'number of columns ({}). Singular '
+                               'values may decrease too slowly.'.format(n))
+        Y[:, j] -= Q[:, :Q_width] @ (np.conj(Q[:, :Q_width]).T @ Y[:, j])
+        q = Y[:, j] / np.linalg.norm(Y[:, j])
+        j = (j + 1) % r
+        Q_width += 1
+        if Q_width > Q.shape[1]:
+            Q_new = np.empty((m, Q.shape[1]*2), dtype=a.dtype)
+            Q_new[:, :Q.shape[1]] = Q
+            Q = Q_new
+        Q[:, Q_width-1] = q
+        # This matrix operation is done on the full matrix and column j-1 is
+        # then overwritten with y to comply with the original algorithm,
+        # even if the orthogonalization of Y would not affect y significantly.
+        Y -= np.outer(q, q.conj() @ Y)
+        w = rand_state.randn(n)
+        A_w = a @ w
+        y = A_w - Q[:, :Q_width] @ (Q[:, :Q_width].T.conj() @ A_w)
+        Y[:, j - 1] = y
+        y_norms[j - 1] = np.linalg.norm(y)
+    return Q[:, :Q_width]
+
+# TODO use similar code for measuring running times
+# def main_loop():
+#     from time import process_time
+#     m = 1000
+#     n = m
+#     p_a = 10
+#     m, n, p_a = 100, 100, 2
+#     a_mat = build_test_matrix(m, n, p=p_a, rand_state=0)
+#     a_op = aslinearoperator(a_mat)
+#     from pomad.utils import FourierMultiplierOp
+#     n = 123
+#     a_op = FourierMultiplierOp(n=n, p=30)
+#     a_mat = a_op @ np.eye(n)
+#     tolerance = 10 ** -3
+#
+#     def adaptive_randomized_range_finder_q2(**nargs):
+#         return adaptive_randomized_range_finder(**nargs, n_cols_Q=2)
+#
+#     def adaptive_randomized_range_finder_q16(**nargs):
+#         return adaptive_randomized_range_finder(**nargs, n_cols_Q=16)
+#
+#     def adaptive_randomized_range_finder_q455(**nargs):
+#         return adaptive_randomized_range_finder(**nargs, n_cols_Q=455)
+#
+#     f_list = [adaptive_randomized_range_finder_naive,
+#               adaptive_randomized_range_finder_nolist,
+#               adaptive_randomized_range_finder_circy,
+#               adaptive_randomized_range_finder_circynorms,
+#               adaptive_randomized_range_finder_rmloop,
+#               adaptive_randomized_range_finder_rmconcQ,
+#               adaptive_randomized_range_finder_mem_alloc,
+#               adaptive_randomized_range_finder_q2,
+#               adaptive_randomized_range_finder_q16,
+#               adaptive_randomized_range_finder_q455,
+#               adaptive_randomized_range_finder,
+#               ]
+#     dt_list = np.empty(len(f_list))
+#     err_list = np.empty(len(f_list))
+#     for i_f, f in enumerate(f_list):
+#         t0 = process_time()
+#         q_mat = f(a=a_op, tolerance=tolerance, r=6, rand_state=0)
+#         dt_list[i_f] = process_time() - t0
+#         err_list[i_f] = np.linalg.norm(a_mat - q_mat @ q_mat.T.conj() @ a_mat,
+#                                        ord=2)
+#         print('t={:.2f}, gain t={:.1%}, err={:.12e}, delta err={:.2e}, '
+#               'Q shape:{} ({})'
+#               .format(dt_list[i_f], dt_list[i_f]/dt_list[0], err_list[i_f],
+#                       abs(err_list[i_f]-err_list[0]), q_mat.shape, f.__name__))
+#
+#
+# if __name__ == '__main__':
+#     main_loop()
diff --git a/pomad/factorization_construction.py b/pomad/factorization_construction.py
new file mode 100755
index 0000000000000000000000000000000000000000..924ea7b51ae341b5de49e0e6dcaa96010da03f43
--- /dev/null
+++ b/pomad/factorization_construction.py
@@ -0,0 +1,144 @@
+# -*- coding: utf-8 -*-
+# ######### COPYRIGHT #########
+# Credits
+# #######
+#
+# Copyright(c) 2019-2020
+# ----------------------
+#
+# * Laboratoire d'Informatique et Systèmes <http://www.lis-lab.fr/>
+# * Université d'Aix-Marseille <http://www.univ-amu.fr/>
+# * Centre National de la Recherche Scientifique <http://www.cnrs.fr/>
+# * Université de Toulon <http://www.univ-tln.fr/>
+#
+# Contributors
+# ------------
+#
+# * Valentin Emiya <firstname.lastname_AT_lis-lab.fr>
+#
+# This package has been created thanks to the joint work with Florent Jaillet
+# and Ronan Hamon on other packages.
+#
+# Description
+# -----------
+#
+# `pomad` is a Python implementation of algorithms from
+# paper *Finding Structure with Randomness: Probabilistic Algorithms for
+# Constructing Approximate Matrix Decompositions*, by N. Halko, P. G.
+# Martinsson and J. A. Tropp, SIAM review, 53 (2), 2011, https://arxiv.org/abs/0909.4061.
+#
+#
+# Version
+# -------
+#
+# * pomad version = 0.1
+#
+# Licence
+# -------
+#
+# 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/>.
+#
+# ######### COPYRIGHT #########
+""" Algorithms for Stage B: Construction of Standard Factorizations
+
+.. moduleauthor:: Valentin Emiya
+"""
+import numpy as np
+from numpy.linalg import cholesky  # TODO check np vs sp implementations
+from scipy.linalg import solve_triangular
+from scipy.sparse.linalg import aslinearoperator, eigs, svds
+
+
+def direct_svd(a_mat, q_mat):
+    """
+    Algo 5.1: direct SVD
+
+    Parameters
+    ----------
+    a_mat : nd-array
+        Positive semidefinite matrix to be decomposed
+    q_mat : nd-array
+        Matrix whose range captures the action of a_mat
+
+    Returns
+    -------
+    u_mat: nd-array
+        Eigenvectors
+    lambda_vec :
+        Eigenvalues
+
+    """
+    b_mat = q_mat.T.conj() @ a_mat
+    u_tilde, lambda_vec, vh_mat = np.linalg.svd(b_mat, full_matrices=False)
+    u_mat = q_mat @ u_tilde
+    return u_mat, lambda_vec, vh_mat
+
+
+def evd_nystrom(a, q_mat):
+    """
+    Algo 5.5: Eigenvalue decomposition via Nyström method
+
+    Parameters
+    ----------
+    a : nd-array or LinearOperator
+        Positive semidefinite matrix or linear operator to be decomposed
+    q_mat : nd-array
+        Matrix whose range captures the action of a_mat
+
+    Returns
+    -------
+    u_mat: nd-array
+        Eigenvectors
+    lambda_vec :
+        Eigenvalues
+    """
+    if isinstance(a, np.ndarray):
+        a = aslinearoperator(a)
+    b1_mat = a(q_mat)
+    b2_mat = q_mat.T.conj() @ b1_mat
+    c = cholesky(b2_mat).T.conj()
+    f_mat = solve_triangular(a=c.T, b=b1_mat.T, lower=True).T
+    u_mat, sigma, _ = np.linalg.svd(f_mat, full_matrices=False)
+    return sigma**2, u_mat
+
+
+# TODO refactor this code into nb
+# if __name__ == '__main__':
+#     from pomad.range_approximation import randomized_range_finder
+#     import matplotlib.pyplot as plt
+#     from pomad.utils import FourierMultiplierOp
+#     print('Test evd with operator')
+#
+#     n = 53
+#     n_l = 41
+#     p = 5
+#     a_op = FourierMultiplierOp(n, p)
+#     q = randomized_range_finder(a_mat=a_op, n_l=n_l)
+#     w_sort = np.sort(a_op.w)[::-1]
+#     plt.semilogy(w_sort, label='w')
+#     print('sing value:', w_sort[n_l])
+#     q_op = aslinearoperator(q)
+#     qh_op = aslinearoperator(q.T.conj())
+#     range_err_op = a_op - q_op @ qh_op @ a_op
+#     print('range approx error:', eigs(range_err_op, k=1)[0][0])
+#     # u51, s51, vh51 = direct_svd(a_mat=a_op, q_mat=q)
+#     # print(eigs(
+#     #     aslinearoperator(SvdErrOp(a_op=a_op, u=u51, s=s51, vh=vh51), k=1)))
+#     s55, u55 = evd_nystrom(a_mat=a_op, q_mat=q)
+#     a_est = u55 @ np.diag(s55) @ u55.T.conj()
+#     print(eigs(a_op - aslinearoperator(a_est))[0][0])
+#
+#     plt.grid()
+#     plt.legend()
+#     e = a_op - aslinearoperator(u55 @ np.diag(s55) @ u55.T.conj())
diff --git a/pomad/range_approximation.py b/pomad/range_approximation.py
new file mode 100755
index 0000000000000000000000000000000000000000..cddba14794a79e5467c7f7ad63e5b5b12a15961a
--- /dev/null
+++ b/pomad/range_approximation.py
@@ -0,0 +1,162 @@
+# -*- coding: utf-8 -*-
+# ######### COPYRIGHT #########
+# Credits
+# #######
+#
+# Copyright(c) 2019-2020
+# ----------------------
+#
+# * Laboratoire d'Informatique et Systèmes <http://www.lis-lab.fr/>
+# * Université d'Aix-Marseille <http://www.univ-amu.fr/>
+# * Centre National de la Recherche Scientifique <http://www.cnrs.fr/>
+# * Université de Toulon <http://www.univ-tln.fr/>
+#
+# Contributors
+# ------------
+#
+# * Valentin Emiya <firstname.lastname_AT_lis-lab.fr>
+#
+# This package has been created thanks to the joint work with Florent Jaillet
+# and Ronan Hamon on other packages.
+#
+# Description
+# -----------
+#
+# `pomad` is a Python implementation of algorithms from
+# paper *Finding Structure with Randomness: Probabilistic Algorithms for
+# Constructing Approximate Matrix Decompositions*, by N. Halko, P. G.
+# Martinsson and J. A. Tropp, SIAM review, 53 (2), 2011, https://arxiv.org/abs/0909.4061.
+#
+#
+# Version
+# -------
+#
+# * pomad version = 0.1
+#
+# Licence
+# -------
+#
+# 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/>.
+#
+# ######### COPYRIGHT #########
+""" Algorithms for Stage A: Randomized Schemes for Approximating the Range
+
+.. moduleauthor:: Valentin Emiya
+"""
+import numpy as np
+
+
+def randomized_range_finder(a, n_l, rand_state=None):
+    """ Algorithm 4.1 : Randomized Range Finder
+
+    Parameters
+    ----------
+    a : nd-array or LinearOperator [m, n]
+        Matrix or linear operator whose range must be computed
+    n_l : int
+        Number of columns in result
+    rand_state : RandomState, int or None
+        If RandomState, random generator.
+        If int or None, random seed used to initialize the pseudo-random
+        number generator.
+
+    Returns
+    -------
+    nd-array [m, n_l]
+        Matrix whose columns are orthonormal and whose range approximates
+        the range of a.
+    """
+    if rand_state is None:
+        rand_state = np.random.RandomState(None)
+    if np.issubdtype(type(rand_state), np.dtype(int).type):
+        rand_state = np.random.RandomState(rand_state)
+    n = a.shape[1]
+    omega = rand_state.randn(n, n_l)
+    y_mat = a @ omega
+    q_mat, _ = np.linalg.qr(y_mat)
+    return q_mat
+
+
+def adaptive_randomized_range_finder(a, tolerance, r=5, proba=None,
+                                     rand_state=None, n_cols_Q=16):
+    """ Algorithm 4.2 : Adaptive Randomized Range Finder
+
+    Parameters
+    ----------
+    a : nd-array or LinearOperator [m, n]
+        Matrix or linear operator whose range must be computed
+    tolerance : float
+        Required tolerance on the range approximation.
+    r : int
+        Number of Gaussian vectors used for a posteriori error estimation.
+        Either `r` or `proba` should be set, not both.
+    proba : float
+        Probability that the error is lower than the tolerance. Either `r` or
+        `proba` should be set, not both.
+    rand_state : RandomState, int or None
+        If RandomState, random generator.
+        If int or None, random seed used to initialize the pseudo-random
+        number generator.
+    n_cols_Q : int
+        Initial number of columns to be allocated in Q.
+
+    Returns
+    -------
+
+    For computational efficiency, this implementation uses circular arrays
+    for `Y` and `y_norms`, memory allocation for `Q` and no internal loop
+    """
+    if proba is None and r is None:
+        raise ValueError('Either proba or n_r should not be None.')
+    if proba is not None and r is not None:
+        raise ValueError('Either proba or n_r should be None.')
+    if proba is not None:
+        r = int(np.ceil(- np.log10((1 - proba) / np.min(a.shape))))
+    if rand_state is None:
+        rand_state = np.random.RandomState(None)
+    if np.issubdtype(type(rand_state), np.dtype(int).type):
+        rand_state = np.random.RandomState(rand_state)
+    m, n = a.shape
+    W = rand_state.randn(n, r)
+    Y = a @ W
+    y_norms = np.linalg.norm(Y, axis=0)
+    j = 0
+    Q_width = 0
+    # order='F' necessary for Q.resize
+    Q = np.empty((m, n_cols_Q), order='F', dtype=a.dtype)
+
+    while np.max(y_norms) > tolerance / (10 * np.sqrt(2 / np.pi)):
+        if Q_width >= n:
+            # TODO to be tested
+            raise RuntimeError('Required tolerance not reached with maximum '
+                               'number of columns ({}). Singular '
+                               'values may decrease too slowly.'.format(n))
+        Y[:, j] -= Q[:, :Q_width] @ (Q[:, :Q_width].T.conj() @ Y[:, j])
+        q = Y[:, j] / np.linalg.norm(Y[:, j])
+        j = (j + 1) % r
+        Q_width += 1
+        if Q_width > Q.shape[1]:
+            Q.resize(m, Q.shape[1]*2)
+        Q[:, Q_width-1] = q
+        # This matrix operation is done on the full matrix and column j-1 is
+        # then overwritten with y to comply with the original algorithm,
+        # even if the orthogonalization of Y would not affect y significantly.
+        Y -= np.outer(q, q.conj() @ Y)
+        w = rand_state.randn(n)
+        A_w = a @ w
+        y = A_w - Q[:, :Q_width] @ (Q[:, :Q_width].T.conj() @ A_w)
+        Y[:, j - 1] = y
+        y_norms[j - 1] = np.linalg.norm(y)
+    return Q[:, :Q_width]
+
diff --git a/pomad/tests/__init__.py b/pomad/tests/__init__.py
new file mode 100755
index 0000000000000000000000000000000000000000..40a96afc6ff09d58a702b76e3f7dd412fe975e26
--- /dev/null
+++ b/pomad/tests/__init__.py
@@ -0,0 +1 @@
+# -*- coding: utf-8 -*-
diff --git a/pomad/tests/test_dev_range_approximation.py b/pomad/tests/test_dev_range_approximation.py
new file mode 100755
index 0000000000000000000000000000000000000000..7a0fc2e453e73048ededa68fccde09983cc54df0
--- /dev/null
+++ b/pomad/tests/test_dev_range_approximation.py
@@ -0,0 +1,119 @@
+# -*- coding: utf-8 -*-
+# ######### COPYRIGHT #########
+# Credits
+# #######
+#
+# Copyright(c) 2019-2020
+# ----------------------
+#
+# * Laboratoire d'Informatique et Systèmes <http://www.lis-lab.fr/>
+# * Université d'Aix-Marseille <http://www.univ-amu.fr/>
+# * Centre National de la Recherche Scientifique <http://www.cnrs.fr/>
+# * Université de Toulon <http://www.univ-tln.fr/>
+#
+# Contributors
+# ------------
+#
+# * Valentin Emiya <firstname.lastname_AT_lis-lab.fr>
+#
+# This package has been created thanks to the joint work with Florent Jaillet
+# and Ronan Hamon on other packages.
+#
+# Description
+# -----------
+#
+# `pomad` is a Python implementation of algorithms from
+# paper *Finding Structure with Randomness: Probabilistic Algorithms for
+# Constructing Approximate Matrix Decompositions*, by N. Halko, P. G.
+# Martinsson and J. A. Tropp, SIAM review, 53 (2), 2011, https://arxiv.org/abs/0909.4061.
+#
+#
+# Version
+# -------
+#
+# * pomad version = 0.1
+#
+# Licence
+# -------
+#
+# 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/>.
+#
+# ######### COPYRIGHT #########
+
+"""Test of the module :module:`pomad._dev_range_approximation`
+
+.. moduleauthor:: Valentin Emiya
+"""
+import unittest
+import numpy as np
+
+from scipy.sparse.linalg import aslinearoperator, svds
+
+from pomad.range_approximation import adaptive_randomized_range_finder
+from pomad.utils import \
+    build_random_psd_matrix, build_test_matrix, FourierMultiplierOp
+from pomad._dev_range_approximation import \
+    adaptive_randomized_range_finder_naive, \
+    adaptive_randomized_range_finder_nolist, \
+    adaptive_randomized_range_finder_circy, \
+    adaptive_randomized_range_finder_circynorms, \
+    adaptive_randomized_range_finder_rmloop, \
+    adaptive_randomized_range_finder_rmconcQ, \
+    adaptive_randomized_range_finder_mem_alloc
+
+
+class TestDevAdaptiveRandomizedRangeFinder(unittest.TestCase):
+    def setUp(self) -> None:
+        m = 79
+        n = 123
+        self.a_dict = {
+            'random matrix':
+                build_test_matrix(m=m, n=n, p=5, rand_state=None),
+            'operator from random matrix':
+                aslinearoperator(build_test_matrix(m=m, n=n, p=5,
+                                                   rand_state=None)),
+            'PSD matrix':
+                build_random_psd_matrix(n=n, rank=m, p=5, rand_state=None),
+            'Fourier multiplier': FourierMultiplierOp(n=n, p=30)
+        }
+        self.f_list = [
+            adaptive_randomized_range_finder_naive,
+            adaptive_randomized_range_finder_nolist,
+            adaptive_randomized_range_finder_circy,
+            adaptive_randomized_range_finder_circynorms,
+            adaptive_randomized_range_finder_rmloop,
+            adaptive_randomized_range_finder_rmconcQ,
+            adaptive_randomized_range_finder_mem_alloc,
+            adaptive_randomized_range_finder,
+        ]
+
+    def test_approximation_error(self):
+        """ Check equation (4.2).
+        """
+        tolerance = 1e-3
+        r = 6
+        for rand_state in (None, 0):
+            for k_a in self.a_dict:
+                err = []
+                a = self.a_dict[k_a]
+                for f in self.f_list:
+                    q = f(a=a, tolerance=tolerance, r=r, rand_state=rand_state)
+
+                    a_op = aslinearoperator(a)
+                    q_op = aslinearoperator(q)
+                    err.append(svds(a_op - q_op @ q_op.H @ a_op, k=1)[1][0])
+                    self.assertLessEqual(err[-1], tolerance, msg='Case ' + k_a)
+                if rand_state is not None: # Use same seed to compare results
+                    np.testing.assert_array_almost_equal(err[:-1], err[1:],
+                                                         err_msg='Case ' + k_a)
diff --git a/pomad/tests/test_factorization_construction.py b/pomad/tests/test_factorization_construction.py
new file mode 100755
index 0000000000000000000000000000000000000000..075bffa1f6f6f16bc97212a0a8eac6705fe1deb2
--- /dev/null
+++ b/pomad/tests/test_factorization_construction.py
@@ -0,0 +1,193 @@
+# -*- coding: utf-8 -*-
+# ######### COPYRIGHT #########
+# Credits
+# #######
+#
+# Copyright(c) 2019-2020
+# ----------------------
+#
+# * Laboratoire d'Informatique et Systèmes <http://www.lis-lab.fr/>
+# * Université d'Aix-Marseille <http://www.univ-amu.fr/>
+# * Centre National de la Recherche Scientifique <http://www.cnrs.fr/>
+# * Université de Toulon <http://www.univ-tln.fr/>
+#
+# Contributors
+# ------------
+#
+# * Valentin Emiya <firstname.lastname_AT_lis-lab.fr>
+#
+# This package has been created thanks to the joint work with Florent Jaillet
+# and Ronan Hamon on other packages.
+#
+# Description
+# -----------
+#
+# `pomad` is a Python implementation of algorithms from
+# paper *Finding Structure with Randomness: Probabilistic Algorithms for
+# Constructing Approximate Matrix Decompositions*, by N. Halko, P. G.
+# Martinsson and J. A. Tropp, SIAM review, 53 (2), 2011, https://arxiv.org/abs/0909.4061.
+#
+#
+# Version
+# -------
+#
+# * pomad version = 0.1
+#
+# Licence
+# -------
+#
+# 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/>.
+#
+# ######### COPYRIGHT #########
+
+"""Test of the module :module:`pomad.factorization_construction`
+
+.. moduleauthor:: Valentin Emiya
+"""
+import unittest
+import numpy as np
+from scipy.sparse.linalg import aslinearoperator, eigs, svds
+
+from pomad.factorization_construction import direct_svd, evd_nystrom
+from pomad.range_approximation import randomized_range_finder
+from pomad.utils import \
+    build_random_psd_matrix, build_test_matrix, FourierMultiplierOp
+
+
+class TestDirectSvd(unittest.TestCase):
+    # TODO add test with complex matrices
+    def test_approximation_error_from_true_range(self):
+        """ Check equation (5.2) using true left singular vectors as matrix Q
+        """
+        m = 47
+        n = 53
+        q_width = 41
+        p = 0
+
+        a = build_test_matrix(m=m, n=n, p=p, rand_state=None)
+        u_true, s_true, vh_true = svds(a, k=q_width)
+
+        # Call to tested code
+        u_est, s_est, vh_est = direct_svd(a_mat=a, q_mat=u_true)
+        
+        eps = np.linalg.norm(a - u_true @ np.diag(s_true) @ vh_true)
+        err = np.linalg.norm(a - u_est @ np.diag(s_est) @ vh_est)
+        try:
+            self.assertLessEqual(err, eps)
+        except AssertionError:
+            np.testing.assert_almost_equal(err - eps, 0)
+
+    def test_approximation_error_from_q_est(self):
+        """ Check equation (5.2) using matrix Q estimated with randomized
+        range finder.
+        """
+        m = 47
+        n = 53
+        q_width = 41
+        p = 0
+
+        a = build_test_matrix(m=m, n=n, p=p, rand_state=None)
+        q = randomized_range_finder(a=a, n_l=q_width)
+
+        # Call to tested code
+        u_est, s_est, vh_est = direct_svd(a_mat=a, q_mat=q)
+
+        eps = np.linalg.norm(a - q @ q.T.conj() @ a)
+        err = np.linalg.norm(a - u_est @ np.diag(s_est) @ vh_est)
+        try:
+            self.assertLessEqual(err, eps)
+        except AssertionError:
+            np.testing.assert_almost_equal(err - eps, 0)
+
+
+class TestEvdNystrom(unittest.TestCase):
+    def setUp(self):
+        n = 53
+        p = 0
+        self.a_dict = {False: build_random_psd_matrix(n=n,
+                                                      p=p,
+                                                      rank=None,
+                                                      rand_state=None,
+                                                      is_complex=False),
+                       True: build_random_psd_matrix(n=n,
+                                                     p=p,
+                                                     rank=None,
+                                                     rand_state=None,
+                                                     is_complex=True)}
+
+    def test_approximation_error_from_true_range(self):
+        """ Check equation (5.2) using true left singular vectors as matrix Q
+        """
+        k = 41
+        for is_complex in (False, True):
+            a = self.a_dict[is_complex]
+            s_true, u_true = eigs(a, k=k)
+            q = u_true
+
+            # Call to tested code
+            s_est, u_est = evd_nystrom(a=a, q_mat=q)
+
+            eps = np.linalg.norm(a - q @ q.T.conj() @ a)
+            err = np.linalg.norm(a - u_est @ np.diag(s_est) @ u_est.T.conj())
+            try:
+                self.assertLessEqual(
+                    err, eps, msg='is_complex={}'.format(is_complex))
+            except AssertionError:
+                np.testing.assert_almost_equal(
+                    err - eps, 0, err_msg='is_complex={}'.format(is_complex))
+
+    def test_approximation_error_from_q_est(self):
+        """ Check equation (5.2) using matrix Q estimated with randomized
+        range finder.
+        """
+        k = 41
+        for is_complex in (False, True):
+            a = self.a_dict[is_complex]
+            q = randomized_range_finder(a=a, n_l=k)
+
+            # Call to tested code
+            s_est, u_est = evd_nystrom(a=a, q_mat=q)
+
+            eps = np.linalg.norm(a - q @ q.T.conj() @ a)
+            err = np.linalg.norm(a - u_est @ np.diag(s_est) @ u_est.T.conj())
+            try:
+                self.assertLessEqual(
+                    err, eps, msg='is_complex={}'.format(is_complex))
+            except AssertionError:
+                np.testing.assert_almost_equal(
+                    err - eps, 0, err_msg='is_complex={}'.format(is_complex))
+
+    def test_approximation_error_operator(self):
+        """ Check equation (5.2) with A as an operator.
+        """
+        n = 53
+        n_l = 41
+        p = 5
+        a_op = FourierMultiplierOp(n, p)
+        q = randomized_range_finder(a=a_op, n_l=n_l)
+
+        # Call to tested code
+        s_est, u_est = evd_nystrom(a=a_op, q_mat=q)
+
+        q_op = aslinearoperator(q)
+        qh_op = aslinearoperator(q.T.conj())
+        eps = eigs(a_op - q_op @ qh_op @ a_op, k=1)[0][0]
+
+        a_op_est = aslinearoperator(u_est @ np.diag(s_est) @ u_est.T.conj())
+        err = eigs(a_op - a_op_est, k=1)[0][0]
+
+        try:
+            self.assertLessEqual(err, eps)
+        except AssertionError:
+            np.testing.assert_almost_equal(err - eps, 0)
diff --git a/pomad/tests/test_range_approximation.py b/pomad/tests/test_range_approximation.py
new file mode 100755
index 0000000000000000000000000000000000000000..70ef880a6924b8b96c2403acbd6dad4045623a07
--- /dev/null
+++ b/pomad/tests/test_range_approximation.py
@@ -0,0 +1,149 @@
+# -*- coding: utf-8 -*-
+# ######### COPYRIGHT #########
+# Credits
+# #######
+#
+# Copyright(c) 2019-2020
+# ----------------------
+#
+# * Laboratoire d'Informatique et Systèmes <http://www.lis-lab.fr/>
+# * Université d'Aix-Marseille <http://www.univ-amu.fr/>
+# * Centre National de la Recherche Scientifique <http://www.cnrs.fr/>
+# * Université de Toulon <http://www.univ-tln.fr/>
+#
+# Contributors
+# ------------
+#
+# * Valentin Emiya <firstname.lastname_AT_lis-lab.fr>
+#
+# This package has been created thanks to the joint work with Florent Jaillet
+# and Ronan Hamon on other packages.
+#
+# Description
+# -----------
+#
+# `pomad` is a Python implementation of algorithms from
+# paper *Finding Structure with Randomness: Probabilistic Algorithms for
+# Constructing Approximate Matrix Decompositions*, by N. Halko, P. G.
+# Martinsson and J. A. Tropp, SIAM review, 53 (2), 2011, https://arxiv.org/abs/0909.4061.
+#
+#
+# Version
+# -------
+#
+# * pomad version = 0.1
+#
+# Licence
+# -------
+#
+# 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/>.
+#
+# ######### COPYRIGHT #########
+
+"""Test of the module :module:`pomad.range_approximation`
+
+.. moduleauthor:: Valentin Emiya
+"""
+import unittest
+import numpy as np
+
+from scipy.sparse.linalg import aslinearoperator, svds
+
+from pomad.range_approximation import \
+    randomized_range_finder, adaptive_randomized_range_finder
+from pomad.utils import \
+    build_random_psd_matrix, build_test_matrix, FourierMultiplierOp
+
+
+# TODO to be implemented: approximation, time, complex values
+class TestRandomizedRangeFinder(unittest.TestCase):
+    def setUp(self) -> None:
+        m = 47
+        n = 53
+        self.a_dict = {
+            'random matrix':
+                build_test_matrix(m=m, n=n, p=0, rand_state=None),
+            'operator from random matrix':
+                aslinearoperator(build_test_matrix(m=m, n=n, p=2,
+                                                   rand_state=None)),
+            'PSD matrix':
+                build_random_psd_matrix(n=n, rank=m, p=0, rand_state=None),
+            'Fourier multiplier': FourierMultiplierOp(n=n)
+        }
+
+    def test_approximation_error(self):
+        """ Check equation (1.9).
+        """
+
+        k = 25
+        p = 5
+        for k_a in self.a_dict:
+            a = self.a_dict[k_a]
+            m, n = a.shape
+
+            # code to be tested
+            q_width = k + p
+            q = randomized_range_finder(a=a, n_l=q_width)
+
+            a_op = aslinearoperator(a)
+            q_op = aslinearoperator(q)
+            _, s_true, _ = svds(a, k=k+1)
+            err = svds(a_op - q_op @ q_op.H @ a_op, k=1)[0][0]
+            err_bound = (1 + 9 * np.sqrt((k + p) * min(m, n))) * s_true[k]
+            self.assertLessEqual(err, err_bound, msg='Case ' + k_a)
+
+
+class TestAdaptiveRandomizedRangeFinder(unittest.TestCase):
+    def setUp(self) -> None:
+        m = 79
+        n = 123
+        self.a_dict = {
+            'random matrix':
+                build_test_matrix(m=m, n=n, p=5, rand_state=None),
+            'operator from random matrix':
+                aslinearoperator(build_test_matrix(m=m, n=n, p=5,
+                                                   rand_state=None)),
+            'PSD matrix':
+                build_random_psd_matrix(n=n, rank=m, p=5, rand_state=None),
+            'Fourier multiplier': FourierMultiplierOp(n=n, p=20)
+        }
+
+    def test_approximation_error(self):
+        """ Check equation (4.2).
+        """
+        # raise NotImplementedError
+        # TODO debug this test, takes too much time (loop in
+        #  adaptive_randomized_range_finder)
+        # k = 11
+        # p = 5
+        tolerance = 1e-4
+        r = 6
+        proba = None
+        rand_state = None
+        n_cols_Q = 4
+        for k_a in self.a_dict:
+            a = self.a_dict[k_a]
+
+            q = adaptive_randomized_range_finder(a=a,
+                                                 tolerance=tolerance,
+                                                 r=r,
+                                                 proba=proba,
+                                                 rand_state=rand_state,
+                                                 n_cols_Q=n_cols_Q
+                                                 )
+
+            a_op = aslinearoperator(a)
+            q_op = aslinearoperator(q)
+            err = svds(a_op - q_op @ q_op.H @ a_op, k=1)[1][0]
+            self.assertLessEqual(err, tolerance, msg='Case ' + k_a)
diff --git a/pomad/tests/test_utils.py b/pomad/tests/test_utils.py
new file mode 100755
index 0000000000000000000000000000000000000000..4c318bde5d81124b942f0fe6f2e86e28f568d189
--- /dev/null
+++ b/pomad/tests/test_utils.py
@@ -0,0 +1,178 @@
+# -*- coding: utf-8 -*-
+# ######### COPYRIGHT #########
+# Credits
+# #######
+#
+# Copyright(c) 2019-2020
+# ----------------------
+#
+# * Laboratoire d'Informatique et Systèmes <http://www.lis-lab.fr/>
+# * Université d'Aix-Marseille <http://www.univ-amu.fr/>
+# * Centre National de la Recherche Scientifique <http://www.cnrs.fr/>
+# * Université de Toulon <http://www.univ-tln.fr/>
+#
+# Contributors
+# ------------
+#
+# * Valentin Emiya <firstname.lastname_AT_lis-lab.fr>
+#
+# This package has been created thanks to the joint work with Florent Jaillet
+# and Ronan Hamon on other packages.
+#
+# Description
+# -----------
+#
+# `pomad` is a Python implementation of algorithms from
+# paper *Finding Structure with Randomness: Probabilistic Algorithms for
+# Constructing Approximate Matrix Decompositions*, by N. Halko, P. G.
+# Martinsson and J. A. Tropp, SIAM review, 53 (2), 2011, https://arxiv.org/abs/0909.4061.
+#
+#
+# Version
+# -------
+#
+# * pomad version = 0.1
+#
+# Licence
+# -------
+#
+# 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/>.
+#
+# ######### COPYRIGHT #########
+
+"""Test of the module :module:`pomad.utils`
+
+.. moduleauthor:: Valentin Emiya
+"""
+import unittest
+
+import numpy as np
+
+from pomad.utils import \
+    build_test_matrix, build_random_psd_matrix, FourierMultiplierOp
+
+
+class TestBuildRandomMatrix(unittest.TestCase):
+    def test_real(self):
+        m, n, p = 13, 17, 1
+        for rand_state in (None, 0):
+            a_mat = build_test_matrix(m=m, n=n, p=p, rand_state=rand_state)
+            np.testing.assert_array_equal(a_mat.shape, (m, n))
+            self.assertFalse(np.iscomplexobj(a_mat))
+
+
+class TestBuildRandomPsdMatrix(unittest.TestCase):
+    def test_real(self):
+        n = 29
+        is_complex = False
+        for p in (0, 1):
+            for rand_state in (None, 0):
+                a_mat = build_random_psd_matrix(n=n,
+                                            p=p,
+                                            rank=None,
+                                            rand_state=rand_state,
+                                            is_complex=is_complex)
+                s_vec, _ = np.linalg.eig(a_mat)
+
+                np.testing.assert_array_equal(a_mat.shape, (n, n))
+                self.assertFalse(np.iscomplexobj(a_mat))
+                np.testing.assert_array_almost_equal(np.imag(s_vec), 0)
+                np.testing.assert_array_equal(np.real(s_vec) > 0, True)
+
+    def test_cpx(self):
+        n = 19
+        is_complex = True
+        for p in (0, 1):
+            a_mat = build_random_psd_matrix(n=n,
+                                            p=p,
+                                            rank=None,
+                                            rand_state=None,
+                                            is_complex=is_complex)
+            s_vec, _ = np.linalg.eig(a_mat)
+
+            np.testing.assert_array_equal(a_mat.shape, (n, n))
+            self.assertTrue(np.iscomplexobj(a_mat))
+            np.testing.assert_array_almost_equal(np.imag(s_vec), 0)
+            np.testing.assert_array_equal(np.real(s_vec) > 0, True)
+
+    def test_rank(self):
+        n = 17
+        r = 11
+        is_complex = True
+        for p in (0, 1):
+            a_mat = build_random_psd_matrix(n=n,
+                                            p=p,
+                                            rank=r,
+                                            rand_state=None,
+                                            is_complex=is_complex)
+            s_vec, _ = np.linalg.eig(a_mat)
+            s_vec.sort()
+            s_vec = s_vec[::-1]
+
+            np.testing.assert_array_equal(a_mat.shape, (n, n))
+            self.assertTrue(np.iscomplexobj(a_mat))
+            np.testing.assert_array_almost_equal(np.imag(s_vec), 0)
+            np.testing.assert_array_equal(np.real(s_vec[:r]) > 0, True)
+            np.testing.assert_array_almost_equal(np.real(s_vec[r:]), 0)
+            np.testing.assert_array_equal(np.isclose(np.real(s_vec[:r]), 0),
+                                          False)
+
+
+class TestFourierMultiplierOp(unittest.TestCase):
+    def test_fourier_mult_op_vec(self):
+        n = 29
+        for p in (0, 1, 5):
+            op = FourierMultiplierOp(n=n, p=p)
+            np.testing.assert_array_equal(op.shape, (n, n))
+
+            w = op.w
+            x = np.random.randn(n)
+            y = np.fft.ifft(w * np.fft.fft(x))
+            np.testing.assert_array_almost_equal(op(x), y)
+            np.testing.assert_array_almost_equal(op @ x, y)
+            # Auto-adjoint
+            np.testing.assert_array_almost_equal(op.H(x), y)
+            np.testing.assert_array_almost_equal(op.H @ x, y)
+
+    def test_fourier_mult_op_mat1(self):
+        n = 17
+        n_samples = 1
+        for p in (0, 1, 7):
+            op = FourierMultiplierOp(n=n, p=p)
+            np.testing.assert_array_equal(op.shape, (n, n))
+
+            w = op.w
+            x_mat = np.random.randn(n, n_samples)
+            y = np.fft.ifft(w[:, None] * np.fft.fft(x_mat, axis=0), axis=0)
+            np.testing.assert_array_almost_equal(op(x_mat), y)
+            np.testing.assert_array_almost_equal(op @ x_mat, y)
+            # Auto-adjoint
+            np.testing.assert_array_almost_equal(op.H(x_mat), y)
+            np.testing.assert_array_almost_equal(op.H @ x_mat, y)
+
+    def test_fourier_mult_op_mat7(self):
+        n = 22
+        n_samples = 7
+        for p in (0, 1, 4):
+            op = FourierMultiplierOp(n=n, p=p)
+            np.testing.assert_array_equal(op.shape, (n, n))
+
+            w = op.w
+            x_mat = np.random.randn(n, n_samples)
+            y = np.fft.ifft(w[:, None] * np.fft.fft(x_mat, axis=0), axis=0)
+            np.testing.assert_array_almost_equal(op(x_mat), y)
+            np.testing.assert_array_almost_equal(op @ x_mat, y)
+            # Auto-adjoint
+            np.testing.assert_array_almost_equal(op.H(x_mat), y)
+            np.testing.assert_array_almost_equal(op.H @ x_mat, y)
diff --git a/pomad/utils.py b/pomad/utils.py
new file mode 100755
index 0000000000000000000000000000000000000000..c7e537d29bcdf4fbc92b5bddb5816476b1852f44
--- /dev/null
+++ b/pomad/utils.py
@@ -0,0 +1,122 @@
+# -*- coding: utf-8 -*-
+# ######### COPYRIGHT #########
+# Credits
+# #######
+#
+# Copyright(c) 2019-2020
+# ----------------------
+#
+# * Laboratoire d'Informatique et Systèmes <http://www.lis-lab.fr/>
+# * Université d'Aix-Marseille <http://www.univ-amu.fr/>
+# * Centre National de la Recherche Scientifique <http://www.cnrs.fr/>
+# * Université de Toulon <http://www.univ-tln.fr/>
+#
+# Contributors
+# ------------
+#
+# * Valentin Emiya <firstname.lastname_AT_lis-lab.fr>
+#
+# This package has been created thanks to the joint work with Florent Jaillet
+# and Ronan Hamon on other packages.
+#
+# Description
+# -----------
+#
+# `pomad` is a Python implementation of algorithms from
+# paper *Finding Structure with Randomness: Probabilistic Algorithms for
+# Constructing Approximate Matrix Decompositions*, by N. Halko, P. G.
+# Martinsson and J. A. Tropp, SIAM review, 53 (2), 2011, https://arxiv.org/abs/0909.4061.
+#
+#
+# Version
+# -------
+#
+# * pomad version = 0.1
+#
+# Licence
+# -------
+#
+# 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/>.
+#
+# ######### COPYRIGHT #########
+"""Utils classes and functions for pomad.
+
+.. moduleauthor:: Valentin Emiya
+"""
+import numpy as np
+from scipy.sparse.linalg import LinearOperator
+
+
+def build_test_matrix(m, n, p=10, rand_state=None, is_complex=False):
+    if rand_state is None:
+        rand_state = np.random.RandomState(None)
+    if np.issubdtype(type(rand_state), np.dtype(int).type):
+        rand_state = np.random.RandomState(rand_state)
+    if is_complex:
+        raise NotImplementedError()
+    A = rand_state.randn(m, n)
+    A = A / np.linalg.norm(A, ord=2)
+    return np.linalg.matrix_power(A @ A.T.conj(), p) @ A
+
+
+def build_random_psd_matrix(n, p=0, rank=None, rand_state=None,
+                            is_complex=False):
+    """
+    Generates a random positive semidefinite matrix
+
+    Parameters
+    ----------
+    n
+    p
+    rank
+    rand_state
+    is_complex
+
+    Returns
+    -------
+
+    """
+    if rank is None:
+        rank = n
+    if rand_state is None:
+        rand_state = np.random.RandomState(None)
+    if np.issubdtype(type(rand_state), np.dtype(int).type):
+        rand_state = np.random.RandomState(rand_state)
+
+    if is_complex:
+        B = rand_state.randn(n, rank) + 1j*rand_state.randn(n, rank)
+    else:
+        B = rand_state.randn(n, rank)
+    A = B @ B.T.conj()
+    A = A / np.linalg.norm(A, ord=2)
+    return np.linalg.matrix_power(A @ A.T.conj(), p) @ A
+
+
+class FourierMultiplierOp(LinearOperator):
+    def __init__(self, n, p=1):
+        self.w = np.random.rand(n) ** p
+        LinearOperator.__init__(self,
+                                shape=(self.w.size, self.w.size),
+                                dtype=np.complex_)
+
+    def _matvec(self, x):
+        if x.ndim == 2:
+            return self._matmat(x)
+        return np.fft.ifft(np.fft.fft(x) * self.w)
+
+    def _matmat(self, x):
+        return np.fft.ifft(np.fft.fft(x, axis=0) * self.w[:, None], axis=0)
+
+    def _adjoint(self):
+        return self
diff --git a/requirements/defaults.txt b/requirements/defaults.txt
new file mode 100755
index 0000000000000000000000000000000000000000..1655484d4aeca0196a4ff03d10c1a61523266bed
--- /dev/null
+++ b/requirements/defaults.txt
@@ -0,0 +1,3 @@
+--index-url https://pypi.python.org/simple/
+
+numpy>=1.13
diff --git a/requirements/dev.txt b/requirements/dev.txt
new file mode 100755
index 0000000000000000000000000000000000000000..b287a430e410a343d9d1e4dbc663b158dbe15fa5
--- /dev/null
+++ b/requirements/dev.txt
@@ -0,0 +1,6 @@
+--index-url https://pypi.python.org/simple/
+
+coverage
+pytest
+pytest-cov
+pytest-randomly
diff --git a/requirements/doc.txt b/requirements/doc.txt
new file mode 100755
index 0000000000000000000000000000000000000000..3f09571e01be27274a79d20a1ca47f9f76be206f
--- /dev/null
+++ b/requirements/doc.txt
@@ -0,0 +1,5 @@
+--index-url https://pypi.python.org/simple/
+
+nbsphinx
+numpydoc
+sphinx
diff --git a/setup.cfg b/setup.cfg
new file mode 100755
index 0000000000000000000000000000000000000000..e793c6e647787d5d3f3d8d375ee260fe38040ac4
--- /dev/null
+++ b/setup.cfg
@@ -0,0 +1,30 @@
+[tool:pytest]
+testpaths = pomad
+addopts = --verbose
+          --cov-report=term-missing
+          --cov-report=html
+          --cov=pomad
+          --doctest-modules
+
+[coverage:run]
+branch = True
+source = pomad
+include = */pomad/*
+omit = */tests/*
+
+[coverage:report]
+exclude_lines =
+    pragma: no cover
+    if self.debug:
+    if settings.DEBUG
+    raise AssertionError
+    raise NotImplementedError
+    if 0:
+    if __name__ == .__main__.:
+    if obj is None: return
+    if verbose > 0:
+    if self.verbose > 0:
+    if verbose > 1:
+    if self.verbose > 1:
+    pass
+    def __str__(self):
diff --git a/setup.py b/setup.py
new file mode 100755
index 0000000000000000000000000000000000000000..cdf66d11d7fd1327626e0f77253fdd0e55017f61
--- /dev/null
+++ b/setup.py
@@ -0,0 +1,144 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+# ######### COPYRIGHT #########
+# Credits
+# #######
+#
+# Copyright(c) 2019-2020
+# ----------------------
+#
+# * Laboratoire d'Informatique et Systèmes <http://www.lis-lab.fr/>
+# * Université d'Aix-Marseille <http://www.univ-amu.fr/>
+# * Centre National de la Recherche Scientifique <http://www.cnrs.fr/>
+# * Université de Toulon <http://www.univ-tln.fr/>
+#
+# Contributors
+# ------------
+#
+# * Valentin Emiya <firstname.lastname_AT_lis-lab.fr>
+#
+# This package has been created thanks to the joint work with Florent Jaillet
+# and Ronan Hamon on other packages.
+#
+# Description
+# -----------
+#
+# `pomad` is a Python implementation of algorithms from
+# paper *Finding Structure with Randomness: Probabilistic Algorithms for
+# Constructing Approximate Matrix Decompositions*, by N. Halko, P. G.
+# Martinsson and J. A. Tropp, SIAM review, 53 (2), 2011, https://arxiv.org/abs/0909.4061.
+#
+#
+# Version
+# -------
+#
+# * pomad version = 0.1
+#
+# Licence
+# -------
+#
+# 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/>.
+#
+# ######### COPYRIGHT #########
+
+import os
+from setuptools import setup, find_packages
+import sys
+
+NAME = 'pomad'
+DESCRIPTION = 'PrObabilistic MAtrix Decompositions from Halko et al., 2011'
+LICENSE = 'GNU General Public License v3 (GPLv3)'
+URL = 'https://gitlab.lis-lab.fr/templates/{}'.format(NAME)
+AUTHOR = 'Valentin Emiya'
+AUTHOR_EMAIL = ('valentin.emiya@lis-lab.fr')
+INSTALL_REQUIRES = ['numpy', 'scipy']
+CLASSIFIERS = [
+    'Development Status :: 5 - Production/Stable',
+    'Intended Audience :: Developers',
+    'Intended Audience :: End Users/Desktop',
+    'Intended Audience :: Science/Research',
+    'Topic :: Scientific/Engineering :: Mathematics',
+    'License :: OSI Approved :: GNU General Public License v3 (GPLv3)',
+    'Natural Language :: English',
+    'Operating System :: MacOS :: MacOS X ',
+    'Operating System :: POSIX :: Linux',
+    'Programming Language :: Python :: 3.5',
+    'Programming Language :: Python :: 3.6']
+PYTHON_REQUIRES = '>=3.5'
+EXTRAS_REQUIRE = {
+    'dev': ['coverage', 'pytest', 'pytest-cov', 'pytest-randomly'],
+    'doc': ['nbsphinx', 'numpydoc', 'sphinx']}
+PROJECT_URLS = {'Bug Reports': URL + '/issues',
+                'Source': URL}
+KEYWORDS = 'template, package'
+
+###############################################################################
+if sys.argv[-1] == 'setup.py':
+    print("To install, run 'python setup.py install'\n")
+
+if sys.version_info[:2] < (3, 5):
+    errmsg = '{} requires Python 3.5 or later ({[0]:d}.{[1]:d} detected).'
+    print(errmsg.format(NAME, sys.version_info[:2]))
+    sys.exit(-1)
+
+
+def get_version():
+    v_text = open('VERSION').read().strip()
+    v_text_formted = '{"' + v_text.replace('\n', '","').replace(':', '":"')
+    v_text_formted += '"}'
+    v_dict = eval(v_text_formted)
+    return v_dict[NAME]
+
+
+def set_version(path, VERSION):
+    filename = os.path.join(path, '__init__.py')
+    buf = ""
+    for line in open(filename, "rb"):
+        if not line.decode("utf8").startswith("__version__ ="):
+            buf += line.decode("utf8")
+    f = open(filename, "wb")
+    f.write(buf.encode("utf8"))
+    f.write(('__version__ = "%s"\n' % VERSION).encode("utf8"))
+
+
+def setup_package():
+    """Setup function"""
+    # set version
+    VERSION = get_version()
+
+    here = os.path.abspath(os.path.dirname(__file__))
+    with open(os.path.join(here, 'README.rst'), encoding='utf-8') as f:
+        long_description = f.read()
+
+    mod_dir = NAME
+    set_version(mod_dir, get_version())
+    setup(name=NAME,
+          version=VERSION,
+          description=DESCRIPTION,
+          long_description=long_description,
+          url=URL,
+          author=AUTHOR,
+          author_email=AUTHOR_EMAIL,
+          license=LICENSE,
+          classifiers=CLASSIFIERS,
+          keywords=KEYWORDS,
+          packages=find_packages(exclude=['doc', 'dev']),
+          install_requires=INSTALL_REQUIRES,
+          python_requires=PYTHON_REQUIRES,
+          extras_require=EXTRAS_REQUIRE,
+          project_urls=PROJECT_URLS)
+
+
+if __name__ == "__main__":
+    setup_package()