Index: /trunk/anuga_core/source/anuga_parallel/pypar_extras/LICENSE
===================================================================
--- /trunk/anuga_core/source/anuga_parallel/pypar_extras/LICENSE	(revision 8514)
+++ /trunk/anuga_core/source/anuga_parallel/pypar_extras/LICENSE	(revision 8514)
@@ -0,0 +1,674 @@
+                    GNU GENERAL PUBLIC LICENSE
+                       Version 3, 29 June 2007
+
+ Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/>
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+                            Preamble
+
+  The GNU General Public License is a free, copyleft license for
+software and other kinds of works.
+
+  The licenses for most software and other practical works are designed
+to take away your freedom to share and change the works.  By contrast,
+the GNU General Public License is intended to guarantee your freedom to
+share and change all versions of a program--to make sure it remains free
+software for all its users.  We, the Free Software Foundation, use the
+GNU General Public License for most of our software; it applies also to
+any other work released this way by its authors.  You can apply it to
+your programs, too.
+
+  When we speak of free software, we are referring to freedom, not
+price.  Our General Public Licenses are designed to make sure that you
+have the freedom to distribute copies of free software (and charge for
+them if you wish), that you receive source code or can get it if you
+want it, that you can change the software or use pieces of it in new
+free programs, and that you know you can do these things.
+
+  To protect your rights, we need to prevent others from denying you
+these rights or asking you to surrender the rights.  Therefore, you have
+certain responsibilities if you distribute copies of the software, or if
+you modify it: responsibilities to respect the freedom of others.
+
+  For example, if you distribute copies of such a program, whether
+gratis or for a fee, you must pass on to the recipients the same
+freedoms that you received.  You must make sure that they, too, receive
+or can get the source code.  And you must show them these terms so they
+know their rights.
+
+  Developers that use the GNU GPL protect your rights with two steps:
+(1) assert copyright on the software, and (2) offer you this License
+giving you legal permission to copy, distribute and/or modify it.
+
+  For the developers' and authors' protection, the GPL clearly explains
+that there is no warranty for this free software.  For both users' and
+authors' sake, the GPL requires that modified versions be marked as
+changed, so that their problems will not be attributed erroneously to
+authors of previous versions.
+
+  Some devices are designed to deny users access to install or run
+modified versions of the software inside them, although the manufacturer
+can do so.  This is fundamentally incompatible with the aim of
+protecting users' freedom to change the software.  The systematic
+pattern of such abuse occurs in the area of products for individuals to
+use, which is precisely where it is most unacceptable.  Therefore, we
+have designed this version of the GPL to prohibit the practice for those
+products.  If such problems arise substantially in other domains, we
+stand ready to extend this provision to those domains in future versions
+of the GPL, as needed to protect the freedom of users.
+
+  Finally, every program is threatened constantly by software patents.
+States should not allow patents to restrict development and use of
+software on general-purpose computers, but in those that do, we wish to
+avoid the special danger that patents applied to a free program could
+make it effectively proprietary.  To prevent this, the GPL assures that
+patents cannot be used to render the program non-free.
+
+  The precise terms and conditions for copying, distribution and
+modification follow.
+
+                       TERMS AND CONDITIONS
+
+  0. Definitions.
+
+  "This License" refers to version 3 of the GNU General Public License.
+
+  "Copyright" also means copyright-like laws that apply to other kinds of
+works, such as semiconductor masks.
+
+  "The Program" refers to any copyrightable work licensed under this
+License.  Each licensee is addressed as "you".  "Licensees" and
+"recipients" may be individuals or organizations.
+
+  To "modify" a work means to copy from or adapt all or part of the work
+in a fashion requiring copyright permission, other than the making of an
+exact copy.  The resulting work is called a "modified version" of the
+earlier work or a work "based on" the earlier work.
+
+  A "covered work" means either the unmodified Program or a work based
+on the Program.
+
+  To "propagate" a work means to do anything with it that, without
+permission, would make you directly or secondarily liable for
+infringement under applicable copyright law, except executing it on a
+computer or modifying a private copy.  Propagation includes copying,
+distribution (with or without modification), making available to the
+public, and in some countries other activities as well.
+
+  To "convey" a work means any kind of propagation that enables other
+parties to make or receive copies.  Mere interaction with a user through
+a computer network, with no transfer of a copy, is not conveying.
+
+  An interactive user interface displays "Appropriate Legal Notices"
+to the extent that it includes a convenient and prominently visible
+feature that (1) displays an appropriate copyright notice, and (2)
+tells the user that there is no warranty for the work (except to the
+extent that warranties are provided), that licensees may convey the
+work under this License, and how to view a copy of this License.  If
+the interface presents a list of user commands or options, such as a
+menu, a prominent item in the list meets this criterion.
+
+  1. Source Code.
+
+  The "source code" for a work means the preferred form of the work
+for making modifications to it.  "Object code" means any non-source
+form of a work.
+
+  A "Standard Interface" means an interface that either is an official
+standard defined by a recognized standards body, or, in the case of
+interfaces specified for a particular programming language, one that
+is widely used among developers working in that language.
+
+  The "System Libraries" of an executable work include anything, other
+than the work as a whole, that (a) is included in the normal form of
+packaging a Major Component, but which is not part of that Major
+Component, and (b) serves only to enable use of the work with that
+Major Component, or to implement a Standard Interface for which an
+implementation is available to the public in source code form.  A
+"Major Component", in this context, means a major essential component
+(kernel, window system, and so on) of the specific operating system
+(if any) on which the executable work runs, or a compiler used to
+produce the work, or an object code interpreter used to run it.
+
+  The "Corresponding Source" for a work in object code form means all
+the source code needed to generate, install, and (for an executable
+work) run the object code and to modify the work, including scripts to
+control those activities.  However, it does not include the work's
+System Libraries, or general-purpose tools or generally available free
+programs which are used unmodified in performing those activities but
+which are not part of the work.  For example, Corresponding Source
+includes interface definition files associated with source files for
+the work, and the source code for shared libraries and dynamically
+linked subprograms that the work is specifically designed to require,
+such as by intimate data communication or control flow between those
+subprograms and other parts of the work.
+
+  The Corresponding Source need not include anything that users
+can regenerate automatically from other parts of the Corresponding
+Source.
+
+  The Corresponding Source for a work in source code form is that
+same work.
+
+  2. Basic Permissions.
+
+  All rights granted under this License are granted for the term of
+copyright on the Program, and are irrevocable provided the stated
+conditions are met.  This License explicitly affirms your unlimited
+permission to run the unmodified Program.  The output from running a
+covered work is covered by this License only if the output, given its
+content, constitutes a covered work.  This License acknowledges your
+rights of fair use or other equivalent, as provided by copyright law.
+
+  You may make, run and propagate covered works that you do not
+convey, without conditions so long as your license otherwise remains
+in force.  You may convey covered works to others for the sole purpose
+of having them make modifications exclusively for you, or provide you
+with facilities for running those works, provided that you comply with
+the terms of this License in conveying all material for which you do
+not control copyright.  Those thus making or running the covered works
+for you must do so exclusively on your behalf, under your direction
+and control, on terms that prohibit them from making any copies of
+your copyrighted material outside their relationship with you.
+
+  Conveying under any other circumstances is permitted solely under
+the conditions stated below.  Sublicensing is not allowed; section 10
+makes it unnecessary.
+
+  3. Protecting Users' Legal Rights From Anti-Circumvention Law.
+
+  No covered work shall be deemed part of an effective technological
+measure under any applicable law fulfilling obligations under article
+11 of the WIPO copyright treaty adopted on 20 December 1996, or
+similar laws prohibiting or restricting circumvention of such
+measures.
+
+  When you convey a covered work, you waive any legal power to forbid
+circumvention of technological measures to the extent such circumvention
+is effected by exercising rights under this License with respect to
+the covered work, and you disclaim any intention to limit operation or
+modification of the work as a means of enforcing, against the work's
+users, your or third parties' legal rights to forbid circumvention of
+technological measures.
+
+  4. Conveying Verbatim Copies.
+
+  You may convey verbatim copies of the Program's source code as you
+receive it, in any medium, provided that you conspicuously and
+appropriately publish on each copy an appropriate copyright notice;
+keep intact all notices stating that this License and any
+non-permissive terms added in accord with section 7 apply to the code;
+keep intact all notices of the absence of any warranty; and give all
+recipients a copy of this License along with the Program.
+
+  You may charge any price or no price for each copy that you convey,
+and you may offer support or warranty protection for a fee.
+
+  5. Conveying Modified Source Versions.
+
+  You may convey a work based on the Program, or the modifications to
+produce it from the Program, in the form of source code under the
+terms of section 4, provided that you also meet all of these conditions:
+
+    a) The work must carry prominent notices stating that you modified
+    it, and giving a relevant date.
+
+    b) The work must carry prominent notices stating that it is
+    released under this License and any conditions added under section
+    7.  This requirement modifies the requirement in section 4 to
+    "keep intact all notices".
+
+    c) You must license the entire work, as a whole, under this
+    License to anyone who comes into possession of a copy.  This
+    License will therefore apply, along with any applicable section 7
+    additional terms, to the whole of the work, and all its parts,
+    regardless of how they are packaged.  This License gives no
+    permission to license the work in any other way, but it does not
+    invalidate such permission if you have separately received it.
+
+    d) If the work has interactive user interfaces, each must display
+    Appropriate Legal Notices; however, if the Program has interactive
+    interfaces that do not display Appropriate Legal Notices, your
+    work need not make them do so.
+
+  A compilation of a covered work with other separate and independent
+works, which are not by their nature extensions of the covered work,
+and which are not combined with it such as to form a larger program,
+in or on a volume of a storage or distribution medium, is called an
+"aggregate" if the compilation and its resulting copyright are not
+used to limit the access or legal rights of the compilation's users
+beyond what the individual works permit.  Inclusion of a covered work
+in an aggregate does not cause this License to apply to the other
+parts of the aggregate.
+
+  6. Conveying Non-Source Forms.
+
+  You may convey a covered work in object code form under the terms
+of sections 4 and 5, provided that you also convey the
+machine-readable Corresponding Source under the terms of this License,
+in one of these ways:
+
+    a) Convey the object code in, or embodied in, a physical product
+    (including a physical distribution medium), accompanied by the
+    Corresponding Source fixed on a durable physical medium
+    customarily used for software interchange.
+
+    b) Convey the object code in, or embodied in, a physical product
+    (including a physical distribution medium), accompanied by a
+    written offer, valid for at least three years and valid for as
+    long as you offer spare parts or customer support for that product
+    model, to give anyone who possesses the object code either (1) a
+    copy of the Corresponding Source for all the software in the
+    product that is covered by this License, on a durable physical
+    medium customarily used for software interchange, for a price no
+    more than your reasonable cost of physically performing this
+    conveying of source, or (2) access to copy the
+    Corresponding Source from a network server at no charge.
+
+    c) Convey individual copies of the object code with a copy of the
+    written offer to provide the Corresponding Source.  This
+    alternative is allowed only occasionally and noncommercially, and
+    only if you received the object code with such an offer, in accord
+    with subsection 6b.
+
+    d) Convey the object code by offering access from a designated
+    place (gratis or for a charge), and offer equivalent access to the
+    Corresponding Source in the same way through the same place at no
+    further charge.  You need not require recipients to copy the
+    Corresponding Source along with the object code.  If the place to
+    copy the object code is a network server, the Corresponding Source
+    may be on a different server (operated by you or a third party)
+    that supports equivalent copying facilities, provided you maintain
+    clear directions next to the object code saying where to find the
+    Corresponding Source.  Regardless of what server hosts the
+    Corresponding Source, you remain obligated to ensure that it is
+    available for as long as needed to satisfy these requirements.
+
+    e) Convey the object code using peer-to-peer transmission, provided
+    you inform other peers where the object code and Corresponding
+    Source of the work are being offered to the general public at no
+    charge under subsection 6d.
+
+  A separable portion of the object code, whose source code is excluded
+from the Corresponding Source as a System Library, need not be
+included in conveying the object code work.
+
+  A "User Product" is either (1) a "consumer product", which means any
+tangible personal property which is normally used for personal, family,
+or household purposes, or (2) anything designed or sold for incorporation
+into a dwelling.  In determining whether a product is a consumer product,
+doubtful cases shall be resolved in favor of coverage.  For a particular
+product received by a particular user, "normally used" refers to a
+typical or common use of that class of product, regardless of the status
+of the particular user or of the way in which the particular user
+actually uses, or expects or is expected to use, the product.  A product
+is a consumer product regardless of whether the product has substantial
+commercial, industrial or non-consumer uses, unless such uses represent
+the only significant mode of use of the product.
+
+  "Installation Information" for a User Product means any methods,
+procedures, authorization keys, or other information required to install
+and execute modified versions of a covered work in that User Product from
+a modified version of its Corresponding Source.  The information must
+suffice to ensure that the continued functioning of the modified object
+code is in no case prevented or interfered with solely because
+modification has been made.
+
+  If you convey an object code work under this section in, or with, or
+specifically for use in, a User Product, and the conveying occurs as
+part of a transaction in which the right of possession and use of the
+User Product is transferred to the recipient in perpetuity or for a
+fixed term (regardless of how the transaction is characterized), the
+Corresponding Source conveyed under this section must be accompanied
+by the Installation Information.  But this requirement does not apply
+if neither you nor any third party retains the ability to install
+modified object code on the User Product (for example, the work has
+been installed in ROM).
+
+  The requirement to provide Installation Information does not include a
+requirement to continue to provide support service, warranty, or updates
+for a work that has been modified or installed by the recipient, or for
+the User Product in which it has been modified or installed.  Access to a
+network may be denied when the modification itself materially and
+adversely affects the operation of the network or violates the rules and
+protocols for communication across the network.
+
+  Corresponding Source conveyed, and Installation Information provided,
+in accord with this section must be in a format that is publicly
+documented (and with an implementation available to the public in
+source code form), and must require no special password or key for
+unpacking, reading or copying.
+
+  7. Additional Terms.
+
+  "Additional permissions" are terms that supplement the terms of this
+License by making exceptions from one or more of its conditions.
+Additional permissions that are applicable to the entire Program shall
+be treated as though they were included in this License, to the extent
+that they are valid under applicable law.  If additional permissions
+apply only to part of the Program, that part may be used separately
+under those permissions, but the entire Program remains governed by
+this License without regard to the additional permissions.
+
+  When you convey a copy of a covered work, you may at your option
+remove any additional permissions from that copy, or from any part of
+it.  (Additional permissions may be written to require their own
+removal in certain cases when you modify the work.)  You may place
+additional permissions on material, added by you to a covered work,
+for which you have or can give appropriate copyright permission.
+
+  Notwithstanding any other provision of this License, for material you
+add to a covered work, you may (if authorized by the copyright holders of
+that material) supplement the terms of this License with terms:
+
+    a) Disclaiming warranty or limiting liability differently from the
+    terms of sections 15 and 16 of this License; or
+
+    b) Requiring preservation of specified reasonable legal notices or
+    author attributions in that material or in the Appropriate Legal
+    Notices displayed by works containing it; or
+
+    c) Prohibiting misrepresentation of the origin of that material, or
+    requiring that modified versions of such material be marked in
+    reasonable ways as different from the original version; or
+
+    d) Limiting the use for publicity purposes of names of licensors or
+    authors of the material; or
+
+    e) Declining to grant rights under trademark law for use of some
+    trade names, trademarks, or service marks; or
+
+    f) Requiring indemnification of licensors and authors of that
+    material by anyone who conveys the material (or modified versions of
+    it) with contractual assumptions of liability to the recipient, for
+    any liability that these contractual assumptions directly impose on
+    those licensors and authors.
+
+  All other non-permissive additional terms are considered "further
+restrictions" within the meaning of section 10.  If the Program as you
+received it, or any part of it, contains a notice stating that it is
+governed by this License along with a term that is a further
+restriction, you may remove that term.  If a license document contains
+a further restriction but permits relicensing or conveying under this
+License, you may add to a covered work material governed by the terms
+of that license document, provided that the further restriction does
+not survive such relicensing or conveying.
+
+  If you add terms to a covered work in accord with this section, you
+must place, in the relevant source files, a statement of the
+additional terms that apply to those files, or a notice indicating
+where to find the applicable terms.
+
+  Additional terms, permissive or non-permissive, may be stated in the
+form of a separately written license, or stated as exceptions;
+the above requirements apply either way.
+
+  8. Termination.
+
+  You may not propagate or modify a covered work except as expressly
+provided under this License.  Any attempt otherwise to propagate or
+modify it is void, and will automatically terminate your rights under
+this License (including any patent licenses granted under the third
+paragraph of section 11).
+
+  However, if you cease all violation of this License, then your
+license from a particular copyright holder is reinstated (a)
+provisionally, unless and until the copyright holder explicitly and
+finally terminates your license, and (b) permanently, if the copyright
+holder fails to notify you of the violation by some reasonable means
+prior to 60 days after the cessation.
+
+  Moreover, your license from a particular copyright holder is
+reinstated permanently if the copyright holder notifies you of the
+violation by some reasonable means, this is the first time you have
+received notice of violation of this License (for any work) from that
+copyright holder, and you cure the violation prior to 30 days after
+your receipt of the notice.
+
+  Termination of your rights under this section does not terminate the
+licenses of parties who have received copies or rights from you under
+this License.  If your rights have been terminated and not permanently
+reinstated, you do not qualify to receive new licenses for the same
+material under section 10.
+
+  9. Acceptance Not Required for Having Copies.
+
+  You are not required to accept this License in order to receive or
+run a copy of the Program.  Ancillary propagation of a covered work
+occurring solely as a consequence of using peer-to-peer transmission
+to receive a copy likewise does not require acceptance.  However,
+nothing other than this License grants you permission to propagate or
+modify any covered work.  These actions infringe copyright if you do
+not accept this License.  Therefore, by modifying or propagating a
+covered work, you indicate your acceptance of this License to do so.
+
+  10. Automatic Licensing of Downstream Recipients.
+
+  Each time you convey a covered work, the recipient automatically
+receives a license from the original licensors, to run, modify and
+propagate that work, subject to this License.  You are not responsible
+for enforcing compliance by third parties with this License.
+
+  An "entity transaction" is a transaction transferring control of an
+organization, or substantially all assets of one, or subdividing an
+organization, or merging organizations.  If propagation of a covered
+work results from an entity transaction, each party to that
+transaction who receives a copy of the work also receives whatever
+licenses to the work the party's predecessor in interest had or could
+give under the previous paragraph, plus a right to possession of the
+Corresponding Source of the work from the predecessor in interest, if
+the predecessor has it or can get it with reasonable efforts.
+
+  You may not impose any further restrictions on the exercise of the
+rights granted or affirmed under this License.  For example, you may
+not impose a license fee, royalty, or other charge for exercise of
+rights granted under this License, and you may not initiate litigation
+(including a cross-claim or counterclaim in a lawsuit) alleging that
+any patent claim is infringed by making, using, selling, offering for
+sale, or importing the Program or any portion of it.
+
+  11. Patents.
+
+  A "contributor" is a copyright holder who authorizes use under this
+License of the Program or a work on which the Program is based.  The
+work thus licensed is called the contributor's "contributor version".
+
+  A contributor's "essential patent claims" are all patent claims
+owned or controlled by the contributor, whether already acquired or
+hereafter acquired, that would be infringed by some manner, permitted
+by this License, of making, using, or selling its contributor version,
+but do not include claims that would be infringed only as a
+consequence of further modification of the contributor version.  For
+purposes of this definition, "control" includes the right to grant
+patent sublicenses in a manner consistent with the requirements of
+this License.
+
+  Each contributor grants you a non-exclusive, worldwide, royalty-free
+patent license under the contributor's essential patent claims, to
+make, use, sell, offer for sale, import and otherwise run, modify and
+propagate the contents of its contributor version.
+
+  In the following three paragraphs, a "patent license" is any express
+agreement or commitment, however denominated, not to enforce a patent
+(such as an express permission to practice a patent or covenant not to
+sue for patent infringement).  To "grant" such a patent license to a
+party means to make such an agreement or commitment not to enforce a
+patent against the party.
+
+  If you convey a covered work, knowingly relying on a patent license,
+and the Corresponding Source of the work is not available for anyone
+to copy, free of charge and under the terms of this License, through a
+publicly available network server or other readily accessible means,
+then you must either (1) cause the Corresponding Source to be so
+available, or (2) arrange to deprive yourself of the benefit of the
+patent license for this particular work, or (3) arrange, in a manner
+consistent with the requirements of this License, to extend the patent
+license to downstream recipients.  "Knowingly relying" means you have
+actual knowledge that, but for the patent license, your conveying the
+covered work in a country, or your recipient's use of the covered work
+in a country, would infringe one or more identifiable patents in that
+country that you have reason to believe are valid.
+
+  If, pursuant to or in connection with a single transaction or
+arrangement, you convey, or propagate by procuring conveyance of, a
+covered work, and grant a patent license to some of the parties
+receiving the covered work authorizing them to use, propagate, modify
+or convey a specific copy of the covered work, then the patent license
+you grant is automatically extended to all recipients of the covered
+work and works based on it.
+
+  A patent license is "discriminatory" if it does not include within
+the scope of its coverage, prohibits the exercise of, or is
+conditioned on the non-exercise of one or more of the rights that are
+specifically granted under this License.  You may not convey a covered
+work if you are a party to an arrangement with a third party that is
+in the business of distributing software, under which you make payment
+to the third party based on the extent of your activity of conveying
+the work, and under which the third party grants, to any of the
+parties who would receive the covered work from you, a discriminatory
+patent license (a) in connection with copies of the covered work
+conveyed by you (or copies made from those copies), or (b) primarily
+for and in connection with specific products or compilations that
+contain the covered work, unless you entered into that arrangement,
+or that patent license was granted, prior to 28 March 2007.
+
+  Nothing in this License shall be construed as excluding or limiting
+any implied license or other defenses to infringement that may
+otherwise be available to you under applicable patent law.
+
+  12. No Surrender of Others' Freedom.
+
+  If conditions are imposed on you (whether by court order, agreement or
+otherwise) that contradict the conditions of this License, they do not
+excuse you from the conditions of this License.  If you cannot convey a
+covered work so as to satisfy simultaneously your obligations under this
+License and any other pertinent obligations, then as a consequence you may
+not convey it at all.  For example, if you agree to terms that obligate you
+to collect a royalty for further conveying from those to whom you convey
+the Program, the only way you could satisfy both those terms and this
+License would be to refrain entirely from conveying the Program.
+
+  13. Use with the GNU Affero General Public License.
+
+  Notwithstanding any other provision of this License, you have
+permission to link or combine any covered work with a work licensed
+under version 3 of the GNU Affero General Public License into a single
+combined work, and to convey the resulting work.  The terms of this
+License will continue to apply to the part which is the covered work,
+but the special requirements of the GNU Affero General Public License,
+section 13, concerning interaction through a network will apply to the
+combination as such.
+
+  14. Revised Versions of this License.
+
+  The Free Software Foundation may publish revised and/or new versions of
+the GNU General Public License from time to time.  Such new versions will
+be similar in spirit to the present version, but may differ in detail to
+address new problems or concerns.
+
+  Each version is given a distinguishing version number.  If the
+Program specifies that a certain numbered version of the GNU General
+Public License "or any later version" applies to it, you have the
+option of following the terms and conditions either of that numbered
+version or of any later version published by the Free Software
+Foundation.  If the Program does not specify a version number of the
+GNU General Public License, you may choose any version ever published
+by the Free Software Foundation.
+
+  If the Program specifies that a proxy can decide which future
+versions of the GNU General Public License can be used, that proxy's
+public statement of acceptance of a version permanently authorizes you
+to choose that version for the Program.
+
+  Later license versions may give you additional or different
+permissions.  However, no additional obligations are imposed on any
+author or copyright holder as a result of your choosing to follow a
+later version.
+
+  15. Disclaimer of Warranty.
+
+  THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY
+APPLICABLE LAW.  EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT
+HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY
+OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO,
+THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+PURPOSE.  THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM
+IS WITH YOU.  SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF
+ALL NECESSARY SERVICING, REPAIR OR CORRECTION.
+
+  16. Limitation of Liability.
+
+  IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS
+THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY
+GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE
+USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF
+DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD
+PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS),
+EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF
+SUCH DAMAGES.
+
+  17. Interpretation of Sections 15 and 16.
+
+  If the disclaimer of warranty and limitation of liability provided
+above cannot be given local legal effect according to their terms,
+reviewing courts shall apply local law that most closely approximates
+an absolute waiver of all civil liability in connection with the
+Program, unless a warranty or assumption of liability accompanies a
+copy of the Program in return for a fee.
+
+                     END OF TERMS AND CONDITIONS
+
+            How to Apply These Terms to Your New Programs
+
+  If you develop a new program, and you want it to be of the greatest
+possible use to the public, the best way to achieve this is to make it
+free software which everyone can redistribute and change under these terms.
+
+  To do so, attach the following notices to the program.  It is safest
+to attach them to the start of each source file to most effectively
+state the exclusion of warranty; and each file should have at least
+the "copyright" line and a pointer to where the full notice is found.
+
+    <one line to give the program's name and a brief idea of what it does.>
+    Copyright (C) <year>  <name of author>
+
+    This program is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+Also add information on how to contact you by electronic and paper mail.
+
+  If the program does terminal interaction, make it output a short
+notice like this when it starts in an interactive mode:
+
+    <program>  Copyright (C) <year>  <name of author>
+    This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
+    This is free software, and you are welcome to redistribute it
+    under certain conditions; type `show c' for details.
+
+The hypothetical commands `show w' and `show c' should show the appropriate
+parts of the General Public License.  Of course, your program's commands
+might be different; for a GUI interface, you would use an "about box".
+
+  You should also get your employer (if you work as a programmer) or school,
+if any, to sign a "copyright disclaimer" for the program, if necessary.
+For more information on this, and how to apply and follow the GNU GPL, see
+<http://www.gnu.org/licenses/>.
+
+  The GNU General Public License does not permit incorporating your program
+into proprietary programs.  If your program is a subroutine library, you
+may consider it more useful to permit linking proprietary applications with
+the library.  If this is what you want to do, use the GNU Lesser General
+Public License instead of this License.  But first, please read
+<http://www.gnu.org/philosophy/why-not-lgpl.html>.
Index: /trunk/anuga_core/source/anuga_parallel/pypar_extras/__init__.py
===================================================================
--- /trunk/anuga_core/source/anuga_parallel/pypar_extras/__init__.py	(revision 8514)
+++ /trunk/anuga_core/source/anuga_parallel/pypar_extras/__init__.py	(revision 8514)
@@ -0,0 +1,40 @@
+# =============================================================================
+# pypar_ext.py - Extension to Parallel Python using MPI
+# Copyright (C) 2012 Stephen G Roberts
+#              (ANU)
+#
+#    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 2 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 (http://www.gnu.org/copyleft/gpl.html)
+#    for more details.
+#
+#    You should have received a copy of the GNU General Public License
+#    along with this program; if not, write to the Free Software
+#    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
+#
+#
+# Contact address: Stephen.Roberts@anu.edu.au
+#
+# Version: See pypar_ext.__version__
+# =============================================================================
+
+"""The python module pypar_ext.py and the C-extension mpiext.c
+implements some extensions to pypar.
+"""
+
+from pypar_ext import *
+from __metadata__ import __version__, __date__, __author__
+
+#Add path of package to PYTHONPATH to allow C-extension to be loaded
+#import sys
+#sys.path += __path__
+
+
+
+
Index: /trunk/anuga_core/source/anuga_parallel/pypar_extras/__metadata__.py
===================================================================
--- /trunk/anuga_core/source/anuga_parallel/pypar_extras/__metadata__.py	(revision 8514)
+++ /trunk/anuga_core/source/anuga_parallel/pypar_extras/__metadata__.py	(revision 8514)
@@ -0,0 +1,6 @@
+# Meta data for pypar
+__version__ = '1.0'
+__date__ = '13 August 2012'
+__author__ = 'Stephen G Roberts'
+
+
Index: /trunk/anuga_core/source/anuga_parallel/pypar_extras/compile.py
===================================================================
--- /trunk/anuga_core/source/anuga_parallel/pypar_extras/compile.py	(revision 8514)
+++ /trunk/anuga_core/source/anuga_parallel/pypar_extras/compile.py	(revision 8514)
@@ -0,0 +1,206 @@
+"""compile.py - compile Python C-extension
+
+   Commandline usage: 
+     python compile.py <filename>
+
+   Usage from within Python:
+     import compile
+     compile.compile(<filename>,..)
+
+   Ole Nielsen, Oct 2001      
+"""   
+
+import numpy
+def compile(FNs=None, CC=None, LD = None, SFLAG = None, verbose = 1):
+  """compile(FNs=None, CC=None, LD = None, SFLAG = None):
+  
+     Compile FN(s) using compiler CC (e.g. mpicc), 
+     Loader LD and shared flag SFLAG.
+     If CC is absent use default compiler dependent on platform
+     if LD is absent CC is used.
+     if SFLAG is absent platform default is used
+     FNs can be either one filename or a list of filenames
+     In the latter case, the first will be used to name so file.
+  """
+  import os, string, sys, types
+  
+  # Input check
+  #
+  assert not FNs is None, 'No filename provided'
+
+  if not type(FNs) == types.ListType:
+    FNs = [FNs]
+
+  # Determine platform and compiler
+  #
+  if sys.platform == 'sunos5':  # Solaris
+    if CC:
+      compiler = CC
+    else:  
+      compiler = 'gcc'
+    if LD:
+      loader = LD
+    else:  
+      loader = compiler
+    if SFLAG:
+      sharedflag = SFLAG
+    else:  
+      sharedflag = 'G'
+      
+  elif sys.platform == 'osf1V5':  # Compaq AlphaServer
+    if CC:
+      compiler = CC
+    else:  
+      compiler = 'cc'
+    if LD:
+      loader = LD
+    else:  
+      loader = compiler
+    if SFLAG:
+      sharedflag = SFLAG
+    else:  
+      sharedflag = 'shared'    
+      
+  elif sys.platform == 'linux2':  # Linux
+    if CC:
+      compiler = CC
+    else:  
+      compiler = 'gcc'
+    if LD:
+      loader = LD
+    else:  
+      loader = compiler
+    if SFLAG:
+      sharedflag = SFLAG
+    else:  
+      sharedflag = 'shared'    
+      
+  elif sys.platform == 'darwin':  # Mac OS X:
+    if CC:
+      compiler = CC
+    else:  
+      compiler = 'cc'
+    if LD:
+      loader = LD
+    else:  
+      loader = compiler
+    if SFLAG:
+      sharedflag = SFLAG
+    else:  
+      sharedflag = 'bundle -flat_namespace -undefined suppress'
+      
+  else:
+    if verbose: print "Unrecognised platform %s - revert to default"\
+                %sys.platform
+    
+    if CC:
+      compiler = CC
+    else:  
+      compiler = 'cc'
+    if LD:
+      loader = LD
+    else:  
+      loader = 'ld'
+    if SFLAG:
+      sharedflag = SFLAG
+    else:  
+      sharedflag = 'G'
+
+            
+  # Find location of include files
+  #
+  python_include = os.path.join(os.path.join(sys.exec_prefix, 'include'),
+                                'python'+sys.version[:3])
+  
+  
+  # Check existence of Python.h
+  #
+  headerfile = python_include + '/Python.h'
+  try:
+    open(headerfile, 'r')
+  except:
+    raise """Did not find Python header file %s.
+    Make sure files for Python C-extensions are installed. 
+    In debian linux, for example, you need to install a
+    package called something like python2.1-dev""" % headerfile
+
+
+  # Get numpy include
+  numpy_include = numpy.get_include()
+  
+  
+  # Check filename(s)
+  #
+  object_files = ''
+  for FN in FNs:        
+    root, ext = os.path.splitext(FN)
+    if ext == '':
+      FN = FN + '.c'
+    elif ext.lower() != '.c':
+      msg = 'Unrecognised extension: ' + FN
+      raise Exception(msg)
+    
+    try:
+      open(FN,'r')
+    except:    
+      msg = 'Could not open: ' + FN
+      raise Exception(msg) 
+
+    if not object_files: root1 = root  # Remember first filename        
+    object_files += root + '.o '  
+  
+  
+    # Compile
+    #
+    s = '%s -c %s -I%s -I%s -o %s.o -Wall' % (compiler, FN,
+                                              python_include, numpy_include,
+                                              root)
+    if os.name == 'posix' and os.uname()[4] == 'x86_64':
+      #Extra flags for 64 bit architectures
+      s += ' -fPIC'
+      
+    
+    if verbose:
+      print s
+    else:
+      s = s + ' 2> /dev/null' # Suppress errors
+  
+    try:
+      os.system(s)
+    except:
+      raise 'Could not compile %s - please try manually' % FN  
+
+  
+  # Make shared library (*.so)
+  s = '%s -%s %s -o %s.so' % (loader, sharedflag, object_files, root1)
+
+  if os.name == 'posix' and os.uname()[4] == 'x86_64':
+      # Extra flags for 64 bit architectures using Portland compilers
+      s += ' -mcmodel=medium'
+  
+  if verbose:
+    print s
+  else:
+    s = s + ' 2> /dev/null' # Suppress warnings
+  
+  try:  
+    err=os.system(s)
+  except:
+    raise 'Could not link %s - please try manually' % root1
+    
+
+
+
+if __name__ == '__main__':
+
+  import sys
+  
+  # Get file to compile
+  #
+  if len(sys.argv) < 2:
+    print 'Usage:'
+    print '  python compile.py <filename>'
+    sys.exit()
+
+  compile(sys.argv[1], verbose=1)
+
Index: /trunk/anuga_core/source/anuga_parallel/pypar_extras/mpiextras.c
===================================================================
--- /trunk/anuga_core/source/anuga_parallel/pypar_extras/mpiextras.c	(revision 8514)
+++ /trunk/anuga_core/source/anuga_parallel/pypar_extras/mpiextras.c	(revision 8514)
@@ -0,0 +1,689 @@
+/************************************************************************/
+/* PyPAR - Parallel Python using MPI                 	                */
+/* Copyright (C) 2001 - 2009 Ole Nielsen                                */
+/*                                                                      */
+/* See enclosed README file for details of installation and use.   	*/
+/*                                                                 	*/   
+/* 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 2 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 (http://www.gnu.org/copyleft/gpl.html)    */
+/* for more details.                                                    */
+/*                                                                      */
+/* You should have received a copy of the GNU General Public License    */
+/* along with this program; if not, write to the Free Software          */
+/*  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307*/
+/*                                                                      */
+/*                                                                 	*/
+/* Contact address: Ole.Moller.Nielsen@gmail.com                        */
+/*                                                                 	*/
+/* version (see __version__ in pypar.py)                                */
+/* date (see __date__ in pypar.py)                                      */
+/************************************************************************/
+
+
+#include "Python.h"
+#include "mpi.h"
+#include "math.h"
+#include "numpy/arrayobject.h"
+
+
+/* Handle MPI constants export (shamelessly stolen from _cursesmodule.c)*/
+#define SetDictInt(string,ch) \
+        PyDict_SetItemString(ModDict, string, PyInt_FromLong((long) (ch)));
+
+/* Remap struct MPI_op to int (easier to handle by python)*/
+#define MAX 1
+#define MIN 2
+#define SUM 3
+#define PROD 4
+#define LAND 5
+#define BAND 6
+#define LOR 7
+#define BOR 8
+#define LXOR 9
+#define BXOR 10
+#define MAXLOC 11
+#define MINLOC 12
+/*#define REPLACE 13 // Not available on all MPI systems */
+
+
+static char errmsg[132];  /*Used to cretae exception messages*/
+
+/* MPI_Bsend() related variables. */
+static void *pt_buf;	/* Pointer to allocated buffer. */
+static int buf_size;	/* Size of buffer to allocate. */ 
+
+int length(PyArrayObject *x) {  
+  /* Compute the total length of contiguous array */
+  /* Necessary for communicating multi dimensional arrays */
+  
+  int i, length;
+  
+  /* Compute total length of contiguous data */
+  length = 1;
+  for (i=0; i<x->nd; i++) {
+    length *= x->dimensions[i];
+  }  
+    
+  return length;
+}  
+
+/**********************************************************/
+/* type_map                                               */
+/*                                                        */
+/**********************************************************/
+MPI_Datatype type_map(PyArrayObject *x, int *count) {  
+  /* Return the MPI Datatype corresponding to                       
+     the Python data type as follows  
+  
+     TYPE    py_type  mpi_type  bytes  symbol
+     ---------------------------------------- 
+     INT       4        6         4      'i'
+     LONG      5        8         8      'l'
+     FLOAT     6       10         4      'f'  
+     DOUBLE    12      11         8      'd'
+  
+     Also return the total number of elements in the array
+  
+     The Python datatype COMPLEX ('F') and COMPLEX_DOUBLE ('D')
+     is treated as a special case to the absence of an 
+     MPI_COMPLEX datatype:
+  
+     Complex arrays are mapped to float or double arrays with real 
+     and imaginary parts alternating and count is updated. */
+  
+  int py_type;
+  MPI_Datatype mpi_type;
+  char err_msg[64];
+
+  *count = length(x);
+  
+  py_type = PyArray_TYPE(x);
+  if (py_type == NPY_DOUBLE) 
+    mpi_type = MPI_DOUBLE;
+  else if (py_type == NPY_INT) 
+    mpi_type = MPI_INT;
+  else if (py_type == NPY_CDOUBLE) {
+    mpi_type = MPI_DOUBLE;
+    (*count) *= 2;
+  } else if (py_type == NPY_FLOAT) 
+    mpi_type = MPI_FLOAT;
+  else if (py_type == NPY_LONG)   
+    mpi_type = MPI_LONG;  
+  else if (py_type == NPY_CFLOAT) {
+    mpi_type = MPI_FLOAT;
+    (*count) *= 2;
+  } else {
+    sprintf(err_msg, 
+	    "Array must be of type int or float. I got py_type == %d", 
+	    py_type);
+    PyErr_SetString(PyExc_ValueError, err_msg);
+    return (MPI_Datatype) NULL;
+  }      
+
+  //printf("Types: py_type=%d, mpi_type=%d\n", py_type, (int) mpi_type);
+  
+  return mpi_type;
+}    
+
+/**********************************************************/
+/* op_map                                                 */
+/*                                                        */
+/**********************************************************/
+MPI_Op op_map(int py_op) {  
+  
+  MPI_Op mpi_op;
+  
+  if (py_op == MAX) 
+    mpi_op = MPI_MAX;
+  else if (py_op == MIN)   
+    mpi_op = MPI_MIN;  
+  else if (py_op == SUM)   
+    mpi_op = MPI_SUM;  
+  else if (py_op == PROD)   
+    mpi_op = MPI_PROD;  
+  else if (py_op == LAND)   
+    mpi_op = MPI_LAND;  
+  else if (py_op == BAND)   
+    mpi_op = MPI_BAND;  
+  else if (py_op == LOR)   
+    mpi_op = MPI_LOR;  
+  else if (py_op == BOR)   
+    mpi_op = MPI_BOR;  
+  else if (py_op == LXOR)   
+    mpi_op = MPI_LXOR;  
+  else if (py_op == BXOR)   
+    mpi_op = MPI_BXOR;  
+  else if (py_op == MAXLOC)
+    mpi_op = MPI_MAXLOC;
+  else if (py_op == MINLOC)   
+    mpi_op = MPI_MINLOC;  
+  /*else if (py_op == REPLACE)*/   
+  /*  mpi_op = MPI_REPLACE; */ 
+  else {
+    PyErr_SetString(PyExc_ValueError, "Operation unknown");
+    return (MPI_Op) NULL;
+  }      
+  
+  return mpi_op;
+}  
+
+
+/**********************************************************/
+/* isend_array (asynchronous)                             */
+/* Send Numpy array of type float, double, int, or long   */
+/*                                                        */
+/**********************************************************/
+static PyObject *isend_array(PyObject *self, PyObject *args) {
+  PyObject *input;
+  PyArrayObject *x;
+  int destination, tag, error, count, myid;
+  MPI_Datatype mpi_type;
+
+  /* process the parameters */
+  if (!PyArg_ParseTuple(args, "Oii", &input, &destination, &tag))
+    return NULL;
+
+  /* Make Numpy array from general sequence type (no cost if already Numpy). */
+  x = (PyArrayObject *)
+    PyArray_ContiguousFromObject(input, NPY_NOTYPE, 0, 0);
+
+  /* Input check and determination of MPI type */
+  mpi_type = type_map(x, &count);
+  if (!mpi_type) return NULL;
+
+  /* call the MPI routine */
+  error = MPI_Send(x->data, count, mpi_type, destination, tag,
+		   MPI_COMM_WORLD);
+  Py_DECREF(x);
+
+  if (error != 0) {
+    MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+    sprintf(errmsg, "Proc %d: MPI_Send failed with error code %d\n",
+	    myid, error);
+    PyErr_SetString(PyExc_RuntimeError, errmsg);
+    return NULL;
+  }
+
+  Py_INCREF(Py_None);
+  return (Py_None);
+}
+
+
+
+
+/*************************************************************/
+/* ireceive_array (asynchronous)                             */
+/* Receive Numpy array of type float, double, int, or long   */
+/*                                                           */
+/*************************************************************/
+static PyObject *ireceive_array(PyObject *self, PyObject *args) {
+  PyArrayObject *x;
+  int source, tag, error, st_length, size, count, myid;
+  MPI_Datatype mpi_type;
+  MPI_Status status;
+
+    
+  if (!PyArg_ParseTuple(args, "Oii", &x, &source, &tag))
+    return NULL;    
+    
+  /* Input check and determination of MPI type */          
+  mpi_type = type_map(x, &count);
+  if (!mpi_type) return NULL;  
+      
+  /* call the MPI routine */
+  error =  MPI_Recv(x->data, count, mpi_type, source, tag,
+		    MPI_COMM_WORLD, &status);
+	 
+  /* Do not DECREF x as it must be returned to Python */
+  	 
+  if (error != 0) {
+    MPI_Comm_rank(MPI_COMM_WORLD, &myid);    
+    sprintf(errmsg, "Proc %d: MPI_Recv failed with error code %d\n", 
+	    myid, error);
+    PyErr_SetString(PyExc_RuntimeError, errmsg);   
+    return NULL;
+  }  
+   
+	 
+  MPI_Get_count(&status, mpi_type, &st_length); 
+  /* status.st_length is not available in all MPI implementations */
+  /* Alternative is: MPI_Get_elements(MPI_Status *, MPI_Datatype, int *); */
+	 
+      
+  /* FIXME: This might not be watertight on all platforms */
+  /* Need C equivalent to itemsize().*/
+  if (mpi_type == MPI_DOUBLE) {
+    size = sizeof(double);  /*8 */
+  } else if (mpi_type == MPI_LONG) {
+    size = sizeof(long); /*8? */
+  } else if (mpi_type == MPI_FLOAT) {
+    size = sizeof(float);
+  } else if (mpi_type == MPI_INT) {    
+    size = sizeof(int);  
+  } else {
+    size = 4;  
+  }
+    
+  return Py_BuildValue("(iiiii)", status.MPI_SOURCE, status.MPI_TAG,
+  status.MPI_ERROR, st_length, size);  
+}
+
+
+
+
+
+/*************************************************************/
+/* allreduce_array                                           */
+/* Allreduce Numpy array of type float, double, int, or long */
+/*                                                           */
+/*************************************************************/
+static PyObject *allreduce_array(PyObject *self, PyObject *args) {
+  PyArrayObject *x;
+  PyArrayObject *d;
+  int op, error, count, count1, myid;
+  MPI_Datatype mpi_type, buffer_type;
+  MPI_Op mpi_op;
+
+  /* process the parameters */
+  if (!PyArg_ParseTuple(args, "OOi", &x, &d, &op)) {
+    PyErr_SetString(PyExc_RuntimeError,
+		    "mpiext.c (allreduce_array): could not parse input");
+    return NULL;
+  }
+
+  /* Input check and determination of MPI type */
+  mpi_type = type_map(x, &count);
+  if (!mpi_type) {
+    PyErr_SetString(PyExc_RuntimeError,
+		    "mpiext.c (allreduce_array): could not determine mpi_type");
+    return NULL;
+  }
+
+  /* This error is caught at the pypar level - so we won't end up here
+     unless mpiext is being used independently */
+  buffer_type = type_map(d, &count1);
+  if (mpi_type != buffer_type) {
+    sprintf(errmsg, "mpiext.c (allreduce_array): Input array and buffer must be of the same type.");
+    PyErr_SetString(PyExc_RuntimeError, errmsg);
+
+    return NULL;
+  }
+
+  if (count != count1) {
+    PyErr_SetString(PyExc_RuntimeError,
+		    "mpiext.c (allreduce_array): Input array and buffer must have same length");
+    return NULL;
+  }
+
+  /* Input check and determination of MPI op */
+  mpi_op = op_map(op);
+  if (!mpi_op) {
+    PyErr_SetString(PyExc_RuntimeError,
+		    "mpiext.c (allreduce_array): could not determine mpi_op");
+    return NULL;
+  }
+
+  if (op == MAXLOC || op == MINLOC) {
+    PyErr_SetString(PyExc_RuntimeError,
+		    "mpiext.c (allreduce_array): MAXLOC and MINLOC are not implemented");
+    return NULL;
+  }
+  else {
+    /* call the MPI routine */
+    error =  MPI_Allreduce(x->data, d->data, count, mpi_type, mpi_op, \
+			MPI_COMM_WORLD);
+  }
+
+  if (error != 0) {
+    MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+    sprintf(errmsg, "Proc %d: MPI_Allreduce failed with error code %d\n",
+	    myid, error);
+    PyErr_SetString(PyExc_RuntimeError, errmsg);
+    return NULL;
+  }
+
+  Py_INCREF(Py_None);
+  return (Py_None);
+}
+
+
+/*************************************************************/
+/* do multiple isends and irecv of Numpy array buffers        */
+/* of type float, double, int, or long                       */
+/*                                                           */
+/*************************************************************/
+static PyObject *sendrecv_array(PyObject *self, PyObject *args) {
+  PyObject *send_bufs;
+  PyArrayObject *send_dest;
+  PyObject *recv_bufs;
+  PyArrayObject *recv_dest;
+
+  PyObject *seq;
+
+  PyArrayObject *x;
+  int op, error, count, count1, myid;
+  MPI_Datatype mpi_type, buffer_type;
+  MPI_Op mpi_op;
+  int i, len, lenx;
+
+  double *xdata;
+
+  /* process the parameters */
+  if (!PyArg_ParseTuple(args, "OOOO", &send_bufs, &send_dest, &recv_bufs, &recv_dest)) {
+    PyErr_SetString(PyExc_RuntimeError,
+		    "mpiextras.c (sendrecv_array): could not parse input");
+    return NULL;
+  }
+
+  seq = PySequence_Fast(send_bufs, "expected a sequence");
+  len = PySequence_Size(send_bufs);
+  printf("len of sequence %d\n",len);
+  for (i = 0; i < len; i++) {
+    x = (PyArrayObject *) PySequence_Fast_GET_ITEM(seq, i);
+    lenx = x->dimensions[0];
+    printf("buf size %d\n",len);
+    xdata = (double *) x->data;
+    printf("x->data[0] %g\n",xdata[0]);
+  }
+  Py_DECREF(seq);
+
+
+/*
+//   Input check and determination of MPI type
+  mpi_type = type_map(x, &count);
+  if (!mpi_type) {
+    PyErr_SetString(PyExc_RuntimeError,
+		    "mpiext.c (allreduce_array): could not determine mpi_type");
+    return NULL;
+  }
+
+
+  //This error is caught at the pypar level - so we won't end up here
+  //  unless mpiext is being used independently
+  buffer_type = type_map(d, &count1);
+  if (mpi_type != buffer_type) {
+    sprintf(errmsg, "mpiext.c (allreduce_array): Input array and buffer must be of the same type.");
+    PyErr_SetString(PyExc_RuntimeError, errmsg);
+
+    return NULL;
+  }
+
+  if (count != count1) {
+    PyErr_SetString(PyExc_RuntimeError,
+		    "mpiext.c (allreduce_array): Input array and buffer must have same length");
+    return NULL;
+  }
+
+  // Input check and determination of MPI op
+  mpi_op = op_map(op);
+  if (!mpi_op) {
+    PyErr_SetString(PyExc_RuntimeError,
+		    "mpiext.c (allreduce_array): could not determine mpi_op");
+    return NULL;
+  }
+
+  if (op == MAXLOC || op == MINLOC) {
+    PyErr_SetString(PyExc_RuntimeError,
+		    "mpiext.c (allreduce_array): MAXLOC and MINLOC are not implemented");
+    return NULL;
+  }
+  else {
+    // call the MPI routine
+    error =  MPI_Allreduce(x->data, d->data, count, mpi_type, mpi_op, \
+			MPI_COMM_WORLD);
+  }
+
+  if (error != 0) {
+    MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+    sprintf(errmsg, "Proc %d: MPI_Allreduce failed with error code %d\n",
+	    myid, error);
+    PyErr_SetString(PyExc_RuntimeError, errmsg);
+    return NULL;
+  }
+*/
+
+
+  Py_INCREF(Py_None);
+  return (Py_None);
+}
+
+
+
+/*************************************************************/
+/* do multiple isends and irecv of Numpy array buffers        */
+/* of type float, double, int, or long                       */
+/*                                                           */
+/*************************************************************/
+static PyObject *send_recv_via_dicts(PyObject *self, PyObject *args) {
+  PyObject *send_dict;
+  PyObject *recv_dict;
+
+  PyObject *seq;
+
+  PyObject *list;
+  PyArrayObject *X;
+  PyArrayObject *Id;
+  int op, error, count, count1, myid;
+  MPI_Datatype mpi_type, buffer_type;
+  MPI_Op mpi_op;
+  int i, k, len, lenx;
+  int num_recv=0;
+  int num_send=0;
+
+  int ierr;
+  MPI_Request requests[20];
+  MPI_Status statuses[20];
+
+  Py_ssize_t pos = 0;
+
+  PyObject *key, *value;
+  long *Iddata;
+  double *xdata;
+
+  /* process the parameters */
+  if (!PyArg_ParseTuple(args, "OO", &send_dict, &recv_dict)) {
+    PyErr_SetString(PyExc_RuntimeError,
+		    "mpiextras.c (sendrecv_array): could not parse input");
+    return NULL;
+  }
+
+  //----------------------------------------------------------------------------
+  // Do the recv first
+  //----------------------------------------------------------------------------
+  num_recv = PyDict_Size(recv_dict);
+  //printf("num_recv = %d\n",num_recv);
+  if (num_recv>10) {
+    PyErr_SetString(PyExc_RuntimeError,
+		    "mpiextras.c; Number of recv communication buffers > 10");
+    return NULL;
+  }
+
+  pos = 0;
+  k = 0;
+  while (PyDict_Next(recv_dict, &pos, &key, &value)) {
+    int i = PyInt_AS_LONG(key);
+    //printf("key %d\n",i);
+
+    Id = (PyArrayObject *) PyList_GetItem(value, 0);
+    X   = (PyArrayObject *) PyList_GetItem(value, 2);
+
+    lenx = X->dimensions[0]*X->dimensions[1];
+
+
+    //printf("buf size %d by 3\n",lenx/3);
+    //xdata = (double *) X->data;
+    //Iddata = (long *) Id->data;
+
+    //printf("k = %d \n",k);
+    ierr = MPI_Irecv(X->data, lenx, MPI_DOUBLE, i, 123, MPI_COMM_WORLD, &requests[k]);
+    //printf("ierr = %d \n",ierr);
+    k++;
+    //for (k=0 ; k < lenx ; k++) printf("X[%d] = %g Id[%d] = %ld \n",k,xdata[k],k,Iddata[k]);
+}
+  //----------------------------------------------------------------------------
+  // Do the sends second
+  //----------------------------------------------------------------------------
+  num_send = PyDict_Size(send_dict);
+  //printf("num_send = %d\n",num_send);
+  if (num_send>10) {
+    PyErr_SetString(PyExc_RuntimeError,
+		    "mpiextras.c; Number of send communication buffers > 10");
+    return NULL;
+  }
+
+  pos = 0;
+  while (PyDict_Next(send_dict, &pos, &key, &value)) {
+    int i = PyInt_AS_LONG(key);
+    //printf("key %d\n",i);
+
+    Id = (PyArrayObject *) PyList_GetItem(value, 0);
+    X   = (PyArrayObject *) PyList_GetItem(value, 2);
+
+    lenx = X->dimensions[0]*X->dimensions[1];
+
+    //printf("buf size %d by 3 \n",lenx/3);
+    //xdata = (double *) X->data;
+    //Iddata = (long *) Id->data;
+    //for (k=0 ; k < lenx ; k++) printf("X[%d] = %g Id[%d] = %ld \n",k,xdata[k],k,Iddata[k]);
+
+    //printf("k = %d \n",k);
+    ierr = MPI_Isend(X->data, lenx, MPI_DOUBLE, i, 123, MPI_COMM_WORLD, &requests[k]);
+    //printf("ierr = %d \n",ierr);
+    
+    k++;
+}
+
+
+  //printf("k = %d\n",k);
+
+
+  //----------------------------------------------------------------------------
+  // Now complete communication. We could put some computation between the 
+  // communication calls above and this call.
+  //----------------------------------------------------------------------------
+  ierr =  MPI_Waitall(k,requests,statuses);
+
+
+
+
+
+/*
+  for (i = 0; i < len; i++) {
+    x = (PyArrayObject *) PySequence_Fast_GET_ITEM(seq, i);
+    lenx = x->dimensions[0];
+    printf("buf size %d\n",len);
+    xdata = (double *) x->data;
+    printf("x->data[0] %g\n",xdata[0]);y
+  }
+  Py_DECREF(seq);
+*/
+
+
+/*
+//   Input check and determination of MPI type
+  mpi_type = type_map(x, &count);
+  if (!mpi_type) {
+    PyErr_SetString(PyExc_RuntimeError,
+		    "mpiext.c (allreduce_array): could not determine mpi_type");
+    return NULL;
+  }
+
+
+  //This error is caught at the pypar level - so we won't end up here
+  //  unless mpiext is being used independently
+  buffer_type = type_map(d, &count1);
+  if (mpi_type != buffer_type) {
+    sprintf(errmsg, "mpiext.c (allreduce_array): Input array and buffer must be of the same type.");
+    PyErr_SetString(PyExc_RuntimeError, errmsg);
+
+    return NULL;
+  }
+
+  if (count != count1) {
+    PyErr_SetString(PyExc_RuntimeError,
+		    "mpiext.c (allreduce_array): Input array and buffer must have same length");
+    return NULL;
+  }
+
+  // Input check and determination of MPI op
+  mpi_op = op_map(op);
+  if (!mpi_op) {
+    PyErr_SetString(PyExc_RuntimeError,
+		    "mpiext.c (allreduce_array): could not determine mpi_op");
+    return NULL;
+  }
+
+  if (op == MAXLOC || op == MINLOC) {
+    PyErr_SetString(PyExc_RuntimeError,
+		    "mpiext.c (allreduce_array): MAXLOC and MINLOC are not implemented");
+    return NULL;
+  }
+  else {
+    // call the MPI routine
+    error =  MPI_Allreduce(x->data, d->data, count, mpi_type, mpi_op, \
+			MPI_COMM_WORLD);
+  }
+
+  if (error != 0) {
+    MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+    sprintf(errmsg, "Proc %d: MPI_Allreduce failed with error code %d\n",
+	    myid, error);
+    PyErr_SetString(PyExc_RuntimeError, errmsg);
+    return NULL;
+  }
+*/
+
+
+  Py_INCREF(Py_None);
+  return (Py_None);
+}
+
+
+ 
+/**********************************/
+/* Method table for python module */
+/**********************************/
+static struct PyMethodDef MethodTable[] = {
+  {"isend_array", isend_array, METH_VARARGS},
+  {"ireceive_array", ireceive_array, METH_VARARGS},
+  {"allreduce_array", allreduce_array, METH_VARARGS},
+  {"sendrecv_array", sendrecv_array, METH_VARARGS},
+  {"send_recv_via_dicts", send_recv_via_dicts, METH_VARARGS},
+  {NULL, NULL}
+};
+
+
+/***************************/
+/* Module initialisation   */
+/***************************/
+void initmpiextras(void){
+  PyObject *m, *ModDict;
+  
+  m = Py_InitModule("mpiextras", MethodTable);
+  
+  /* to handle MPI symbolic constants*/
+  ModDict = PyModule_GetDict(m); 
+  SetDictInt("MPI_ANY_TAG", MPI_ANY_TAG);
+  SetDictInt("MPI_TAG_UB", MPI_TAG_UB);  
+  SetDictInt("MPI_ANY_SOURCE", MPI_ANY_SOURCE);
+  SetDictInt("MAX", MAX);
+  SetDictInt("MIN", MIN);
+  SetDictInt("SUM", SUM);
+  SetDictInt("PROD", PROD);
+  SetDictInt("LAND", LAND);
+  SetDictInt("BAND", BAND);
+  SetDictInt("LOR", LOR);
+  SetDictInt("BOR", BOR);
+  SetDictInt("LXOR", LXOR);
+  SetDictInt("BXOR", BXOR);
+
+   
+  import_array();     /* Necessary for handling of numpy structures  */
+}
Index: /trunk/anuga_core/source/anuga_parallel/pypar_extras/pypar_ext.py
===================================================================
--- /trunk/anuga_core/source/anuga_parallel/pypar_extras/pypar_ext.py	(revision 8514)
+++ /trunk/anuga_core/source/anuga_parallel/pypar_extras/pypar_ext.py	(revision 8514)
@@ -0,0 +1,528 @@
+# =============================================================================
+# pypar_ext.py - Extension to Parallel Python using MPI
+# Copyright (C) 2012 Stephen G Roberts
+#              (ANU)
+#
+#    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 2 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 (http://www.gnu.org/copyleft/gpl.html)
+#    for more details.
+#
+#    You should have received a copy of the GNU General Public License
+#    along with this program; if not, write to the Free Software
+#    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
+#
+#
+# Contact address: Stephen.Roberts@anu.edu.au
+#
+# Version: See pypar_ext.__version__
+# =============================================================================
+
+"""Module pypar_ext.py -Extension to pypar
+
+Need to import pypar first to init mpi.
+
+Public functions:
+
+
+
+isend() -- Asyncronous send (arrays)
+receive() --  Asyncronous receive (arrays)
+allreduce() -- wrapper for MPI_Allreduce (array)
+
+
+See doc strings of individual functions for detailed documentation.
+"""
+
+from numpy import zeros, reshape, product
+from __metadata__ import __version__, __date__, __author__
+
+
+# Constants
+#
+max_tag = 32767      # Max tag value (MPI_TAG_UB didn't work and returned 0)
+control_tag = 13849  # Reserved tag used to identify control information
+default_tag = 1      # Tag used as default if not specified
+
+control_sep = ':'          # Separator for fields in control info (NOT ',')
+control_data_max_size = 64 # Maximal size of string holding control data
+
+
+#---------------------------------------------------------------------------
+# Communication functions
+#--------------------------------------------------------------------------
+
+def isend(x, destination, use_buffer=False, vanilla=False,
+         tag=default_tag, bypass=False):
+    """Wrapper for easy MPI send.
+       Send x to destination.
+
+       Automatically determine appropriate protocol
+       and call corresponding send function.
+       Also passes type and size information on as preceding message to
+       simplify the receive call.
+
+       The variable x can be any (picklable) type, but
+       numpy variables and text strings will most efficient.
+       Setting vanilla = 1 forces vanilla mode for any type.
+
+       If use_buffer is True, workspace x will be used for the return value.
+       In this case the corresponding receive call must specify a buffer.
+       Otherwise a new workspace will be created by receive.
+
+       If bypass is True, all admin and error checks
+       get bypassed to reduce the latency. Should only
+       be used for sending numpy arrays and should be matched
+       with a bypass in the corresponding receive command.
+
+    """
+    import types, string
+
+    if bypass:
+        isend_array(x, destination, tag)
+        return
+
+    # Input check
+    errmsg = 'Destination id (%s) must be an integer.' % destination
+    assert type(destination) == types.IntType, errmsg
+
+    errmsg = 'Tag %d is reserved by pypar - please use another.' % control_tag
+    assert tag != control_tag, errmsg
+
+    # Create metadata about object to be sent
+    control_info, x = create_control_info(x, vanilla, return_object=True)
+    protocol = control_info[0]
+
+
+    # Possibly transmit control data
+    if use_buffer is False:
+        send_control_info(control_info, destination)
+
+
+    # Transmit payload data
+    if protocol == 'array':
+        send_array(x, destination, tag)
+    elif protocol in ['string', 'vanilla']:
+        send_string(x, destination, tag)
+    else:
+        raise 'Unknown protocol: %s' %protocol
+
+
+def ireceive(source, buffer=None, vanilla=False, tag=default_tag,
+            return_status=False, bypass=False):
+    """receive - blocking MPI receive
+
+       Receive data from source.
+
+       Optional parameters:
+         buffer: Use specified buffer for received data (faster). Default None.
+         vanilla: Specify to enforce vanilla protocol for any type.
+                  Default False
+
+         tag: Only received messages tagged as specified. Default default_tag
+         return_status: Return Status object along with result. Default False.
+
+       If no buffer is specified, receive will try to receive a
+       preceding message containing protocol, type, size and shape and
+       then create a suitable buffer.
+
+       If buffer is specified the corresponding send must specify
+       use_buffer = True.
+       The variable buffer can be any (picklable) type, but
+       numpy variables and text strings will most efficient.
+
+       Appropriate protocol will be automatically determined
+       and corresponding receive function called.
+
+
+       If bypass is True, all admin and error checks
+       get bypassed to reduce the latency. Should only
+       be used for receiving numpy arrays and should
+       be matched with a bypass in the corresponding send command.
+       Also buffer must be specified.
+    """
+
+    if bypass:
+        # errmsg = 'bypass mode must be used with specified buffer'
+        # assert buffer is not None, msg
+        stat = ireceive_array(buffer, source, tag)
+    else:
+
+        import types
+
+        # Input check
+        errmsg = 'Source id (%s) must be an integer.' %source
+        assert type(source) == types.IntType, errmsg
+
+        errmsg = 'Tag %d is reserved by pypar - please use another.'\
+            % control_tag
+        assert tag != control_tag, errmsg
+
+
+        # Either receive or create metadata about object to receive
+        if buffer is None:
+            control_info, source = receive_control_info(source,
+                                                        return_source=True)
+            protocol, typecode, size, shape = control_info
+        else:
+            protocol, typecode, size, shape = create_control_info(buffer,
+                                                                  vanilla)
+
+
+        # Receive payload data
+        if protocol == 'array':
+            if buffer is None:
+                buffer = zeros(size,typecode)
+                buffer = reshape(buffer, shape)
+
+            stat = receive_array(buffer, source, tag)
+
+        elif protocol == 'string':
+            if buffer is None:
+                buffer = ' '*size
+
+            stat = receive_string(buffer, source, tag)
+
+        elif protocol == 'vanilla':
+            from cPickle import dumps, loads, UnpicklingError
+            if buffer is None:
+                s = ' '*size
+            else:
+                s = dumps(buffer, protocol=2)
+                s = s + ' '*int(0.1*len(s)) #safety
+
+            stat = receive_string(s, source, tag)
+            try:
+                buffer = loads(s)   #Replace buffer with received result
+            except UnpicklingError, err:
+                raise UnpicklingError(str(err) + " - '%s'" % s)
+        else:
+            raise 'Unknown protocol: %s' %protocol
+
+    # Return received data and possibly the status object
+    if return_status:
+        return buffer, Status(stat)
+    else:
+        return buffer
+
+
+
+
+def allreduce(x, op, buffer=None, vanilla=0, bypass=False):
+    """Allreduce elements in x to buffer (of the same size as x)
+       applying operation op elementwise.
+
+       If bypass is True, all admin and error checks
+       get bypassed to reduce the latency.
+       The buffer must be specified explicitly in this case.
+    """
+
+    if bypass:
+        allreduce_array(x, buffer, op)
+        return
+
+
+    import types
+    from mpiext import size
+    numproc = size()         # Needed to determine buffer size
+
+
+
+    # Create metadata about object
+    protocol, typecode, size, shape = create_control_info(x)
+
+    # Allreduce
+    if protocol == 'array':
+        if buffer is None:
+            buffer = zeros(size*numproc, typecode)
+
+            # Modify shape along axis=0 to match size
+            shape = list(shape)
+            shape[0] *= numproc
+            buffer = reshape(buffer, shape)
+
+
+        msg = 'Data array and buffer must have same type '
+        msg = 'in allreduce. I got types "%s" and "%s"' % (x.dtype.char,
+                                                        buffer.dtype.char)
+        assert x.dtype.char == buffer.dtype.char, msg
+        allreduce_array(x, buffer, op)
+
+
+    elif (protocol == 'vanilla' or protocol == 'string'):
+        raise 'Protocol: %s unsupported for allreduce' % protocol
+    else:
+        raise 'Unknown protocol: %s' % protocol
+
+    return buffer
+
+
+
+
+#---------------------------------------------------------
+# INTERNAL FUNCTIONS
+#---------------------------------------------------------
+
+class Status:
+    """ MPI Status block returned by receive if
+        specified with parameter return_status=True
+    """
+
+    def __init__(self, status_tuple):
+        self.source = status_tuple[0]  # Id of sender
+        self.tag = status_tuple[1]     # Tag of received message
+        self.error = status_tuple[2]   # MPI Error code
+        self.length = status_tuple[3]  # Number of elements transmitted
+        self.size = status_tuple[4]    # Size of one element
+
+    def __repr__(self):
+        return 'Pypar Status Object:\n  source=%d\n  tag=%d\n  error=%d\n  length=%d\n  size=%d\n' %(self.source, self.tag, self.error, self.length, self.size)
+
+    def bytes(self):
+        """Number of bytes transmitted (excl control info)
+        """
+        return self.length * self.size
+
+
+
+def create_control_info(x, vanilla=0, return_object=False):
+    """Determine which protocol to use for communication:
+       (numpy) arrays, strings, or vanilla based x's type.
+
+       There are three protocols:
+       'array':   numpy arrays of type 'i', 'l', 'f', 'd', 'F' or 'D' can be
+                  communicated  with mpiext.send_array and mpiext.receive_array.
+       'string':  Text strings can be communicated with mpiext.send_string and
+                  mpiext.receive_string.
+       'vanilla': All other types can be communicated as string representations
+                  provided that the objects
+                  can be serialised using pickle (or cPickle).
+                  The latter mode is less efficient than the
+                  first two but it can handle general structures.
+
+       Rules:
+       If keyword argument vanilla == 1, vanilla is chosen regardless of
+       x's type.
+       Otherwise if x is a string, the string protocol is chosen
+       If x is an array, the 'array' protocol is chosen provided that x has one
+       of the admissible typecodes.
+
+       The optional argument return_object asks to return object as well.
+       This is useful in case it gets modified as in the case of general structures
+       using the vanilla protocol.
+    """
+
+    import types
+
+    # Default values
+    protocol = 'vanilla'
+    typecode = ' '
+    size = 0
+    shape = ()
+
+    # Determine protocol in case
+    if not vanilla:
+        if type(x) == types.StringType:
+            protocol = 'string'
+            typecode = 'c'
+            size = len(x)
+        elif type(x).__name__ == 'ndarray': # numpy isn't imported yet
+            try:
+                import numpy
+            except:
+                print "WARNING (pypar.py): numpy module could not be imported,",
+                print "reverting to vanilla mode"
+                protocol = 'vanilla'
+            else:
+                typecode = x.dtype.char
+                if typecode in ['i', 'l', 'f', 'd', 'F', 'D']:
+                    protocol = 'array'
+                    shape = x.shape
+                    size = product(shape)
+                else:
+                    print "WARNING (pypar.py): numpy object type %s is not supported."\
+                          %(x.dtype.char)
+                    print "Only types 'i', 'l', 'f', 'd', 'F', 'D' are supported,",
+                    print "Reverting to vanilla mode."
+                    protocol = 'vanilla'
+
+    # Pickle general structures using the vanilla protocol
+    if protocol == 'vanilla':
+        from cPickle import dumps
+        x = dumps(x, protocol=2)
+        size = len(x) # Let count be length of pickled object
+
+    # Return
+    if return_object:
+        return [protocol, typecode, size, shape], x
+    else:
+        return [protocol, typecode, size, shape]
+
+
+
+#----------------------------------------------
+
+def send_control_info(control_info, destination):
+    """Send control info to destination
+    """
+    import string
+
+    # Convert to strings
+    control_info = [str(c) for c in control_info]
+
+    control_msg = string.join(control_info,control_sep)
+    if len(control_msg) > control_data_max_size:
+        errmsg = 'Length of control_info exceeds specified maximium (%d)'\
+                 %control_data_max_size
+        errmsg += ' - Please increase it (in pypar.py)'
+        raise errmsg
+
+    send_string(control_msg, destination, control_tag)
+
+
+def receive_control_info(source, return_source=False):
+    """Receive control info from source
+
+    The optional argument (due to Jim Bosch) also returns the actual source node
+    which can be used to require that the data message come from the same node.
+    """
+
+    # FIXME (Ole): Perhaps we should include actual source in the control info?
+
+    import string
+
+    msg = ' '*control_data_max_size
+
+    stat = receive_string(msg, source, control_tag)
+    # No need to create status object here - it is reserved
+    # for payload communications only
+
+    msg = msg[:stat[3]] # Trim buffer to actual received length (needed?)
+
+    control_info = msg.split(control_sep)
+
+    assert len(control_info) == 4, 'len(control_info) = %d' %len(control_info)
+    control_info[2] = eval(control_info[2]) # Convert back to int
+    control_info[3] = eval(control_info[3]) # Convert back to tuple
+
+    if return_source:
+        return control_info, int(stat[0])
+    else:
+        return control_info
+
+
+#----------------------------------------------------------------------------
+# Initialise module
+#----------------------------------------------------------------------------
+
+
+# Take care of situation where module is part of package
+import sys, os, string, os.path
+dirname = os.path.dirname(string.replace(__name__,'.',os.sep)).strip()
+
+if not dirname:
+    dirname = '.'
+
+if dirname[-1] != os.sep:
+    dirname += os.sep
+
+
+
+# Import MPI extension
+#
+# Verify existence of mpiext.so.
+
+try:
+    import mpiextras
+except:
+    errmsg = 'ERROR: C extension mpiextras could not be imported.\n'
+    errmsg += 'Please compile mpiextras.c manually.\n'
+    #raise Exception, errmsg
+    error = 1
+    print errmsg
+else:
+
+    # Determine if MPI program is allowed to run sequentially on current platform
+    # Attempting to check this automatically may case some systems to hang.
+
+    if sys.platform in ['linux2', 'sunos5', 'win32', 'darwin']:
+        # Linux (LAM,MPICH) or Sun (MPICH)
+        error = 0  #Sequential execution of MPI is allowed
+    else:
+        # Platform: Alpha 'osf1V5'
+        cmdstring = '"import mpiext, sys; mpiext.init(sys.argv); mpiext.finalize()"'
+        #s = 'cd %s; python -c %s' %(dirname, cmdstring)
+        s = 'python -c %s >/dev/null 2>/dev/null' %cmdstring
+        error = os.system(s)
+
+        # The check is performed in a separate shell.
+        # Reason: The Alpha server, LAM/Linux or the Sun cannot recover from a
+        # try:
+        #   mpiext.init(sys.argv)
+
+        # However, on LAM/Linux, this test causes system to hang.
+        # Verified (OMN 12/12/2)
+        # If lamboot is started, the system, will hang when init is called
+        # again further down in this file.
+        # If lamboot is not loaded error will be nozero as it should.
+        # I don't know how to deal with this
+        #
+        #Comparisons of two strategies using LAM
+        #
+        # Strategy 1: Assume seq execution is OK (i.e. set error = 0)
+        # Strategy 2: Try to test if mpi can be initialised (in a separate shell)
+        #
+        #
+        # Strategy 1 (currently used)
+        #                    | Lam booted  | Lam not booted
+        #-----------------------------------------------------
+        #
+        # Sequential exec    |  OK         | Not OK
+        # Parallel exec      |  OK         | Not OK
+        #
+        #
+        # Strategy 2
+        #                    | Lam booted  | Lam not booted
+        #-----------------------------------------------------
+        #
+        # Sequential exec    |  Hangs      | Not OK
+        # Parallel exec      |  Hangs      | OK
+        #
+
+
+
+# Initialise MPI
+#
+# Attempt to initialise mpiext.so
+# If this fails, define a rudimentary interface suitable for
+# sequential execution.
+
+if error:
+    print "WARNING: MPI library could not be initialised"
+
+else:
+    from mpiextras import \
+         isend_array, \
+         ireceive_array, \
+         allreduce_array
+
+    # Work around bug in OpenMPI (December 2009):
+    # https://bugs.launchpad.net/ubuntu/+source/petsc4py/+bug/232036
+
+    from ctypes import *
+    CDLL('libmpi.so', RTLD_GLOBAL)
+    # End work around
+
+
+
+
+
+
+
+
Index: /trunk/anuga_core/source/anuga_parallel/pypar_extras/setup.py
===================================================================
--- /trunk/anuga_core/source/anuga_parallel/pypar_extras/setup.py	(revision 8514)
+++ /trunk/anuga_core/source/anuga_parallel/pypar_extras/setup.py	(revision 8514)
@@ -0,0 +1,200 @@
+# Examples of setup
+#
+# python setup.py install
+# python setup.py install --prefix=/opt/python-2.3
+# python setup.py install --home=~
+#
+# Some times you'll have to add a file called pypar.pth
+# containing the word pypar to site-packages
+#
+# See http://docs.python.org/dist/pure-pkg.html for more on distutils
+
+# FIXME: Now mpiext.c and pypar.py are assumed to be in this directory.
+# Maybe, we should put them in the default package directory, pypar.
+# The repository structure would then be
+#
+# pypar
+#     demos
+#     documentation
+#     source
+#          pypar
+
+from distutils.core import setup, Extension
+
+# FIXME (Ole): This works, but I don't know how to use it
+# Generate Python EGG if possible.
+#try:
+#   from setuptools import setup, Extension
+#except ImportError:
+#   pass
+
+import distutils.sysconfig
+import distutils.debug
+import os, sys
+import popen2 # FIXME: Replace with subprocess - http://docs.python.org/library/subprocess.html#replacing-older-functions-with-the-subprocess-module
+import string
+import tempfile
+import numpy
+import shlex
+from __metadata__ import __version__, __date__, __author__
+
+
+def setup_compiler():
+    distutils.sysconfig.get_config_vars()
+    config_vars = distutils.sysconfig._config_vars
+
+    if sys.platform == 'sunos5':
+        config_vars['LDSHARED'] = 'gcc -G'
+        config_vars['CCSHARED'] = ''
+
+
+def uniq_arr(arr):
+    """Remove repeated values from an array and return new array."""
+    ret = []
+    for i in arr:
+        if i not in ret:
+            ret.append(i)
+    return ret
+
+
+def _run_command(cmd):
+    out_file, in_file, err_file = popen2.popen3(cmd)
+    output = out_file.read() + err_file.read()
+    out_file.close()
+    in_file.close()
+    err_file.close()
+    # need this hack to get the exit status
+    print 'running ' + cmd
+    out_file = os.popen(cmd)
+    if out_file.close():
+        # close returns exit status of command.
+        return ''
+    else:
+        # no errors, out_file.close() returns None.
+        return output
+
+
+def _get_mpi_cmd():
+    """Returns the output of the command used to compile using
+    mpicc."""
+    # LAM/OPENMPI/MPICH2
+    output = _run_command('mpicc -show')
+    if output:
+        return output
+
+    # FIXME: If appears that MPICH actually needs this these hacks.
+
+    # MPICH
+    # works with MPICH version 1.2.1 (on Debian)
+    output = _run_command('mpicc -compile_info -link_info')
+    if output:
+        return output
+
+    # Old version of MPICH needs this hack.
+    tmp_base = tempfile.mktemp()
+    tmp_c = tmp_base + ".c"
+    tmp_o = tmp_base + ".o"
+    tmp_file = open(tmp_c, "w")
+    tmp_file.write('#include "mpi.h"\nint main(){return 0;}\n')
+    tmp_file.close()
+    output = _run_command("mpicc -show;"\
+                          "mpicc -echo -c %s -o %s"%(tmp_c, tmp_o))
+    os.remove(tmp_c)
+    if os.path.exists(tmp_o):
+        os.remove(tmp_o)
+    if output:
+        return output
+    else:
+        return ''
+
+
+def get_mpi_flags():
+    output = _get_mpi_cmd()
+    print output
+    if not output:
+        if sys.platform=='win32': # From Simon Frost
+            #this didn't work on my machine (Vladimir Lazunin on April 7, 2009)
+            #output = "gcc -L$MPICH_DIR\SDK.gcc\lib -lmpich -I$MPICH_DIR\SDK.gcc\include"
+
+            #"MPICH_DIR" must be set manually in environment variables
+            mpi_dir = os.getenv("MPICH_DIR")
+            if mpi_dir == None:
+            	print 'MPICH_DIR environment variable must be set'
+            	exit()
+
+            #for MPICH2
+            sdk_prefix = mpi_dir
+            lib_name = 'mpi'
+
+            #for MPICH1
+            if os.path.exists(sdk_prefix + '\\SDK'):
+                sdk_prefix += '\\SDK'
+                lib_name = 'mpich'
+            output = 'gcc -L"%(sdk_prefix)s\lib" -l"%(lib_name)s" -I"%(sdk_prefix)s\include"' % {'sdk_prefix' : sdk_prefix, 'lib_name' : lib_name}
+        else:
+            output = 'cc -L/usr/opt/mpi -lmpi -lelan'
+
+
+    # Now get the include, library dirs and the libs to link with.
+    #flags = string.split(output)
+    flags = shlex.split(output)
+    flags = uniq_arr(flags) # Remove repeated values.
+    inc_dirs = []
+    lib_dirs = []
+    libs = []
+    def_macros = []
+    undef_macros = []
+    for f in flags:
+        if f[:2] == '-I':
+            inc_dirs.append(f[2:])
+        elif f[:2] == '-L':
+            lib_dirs.append(f[2:])
+        elif f[:2] == '-l' and f[-1] != "'": # Patched by Michael McKerns July 2009
+            libs.append(f[2:])
+        elif f[:2] == '-U':
+            undef_macros.append(f[2:])
+        elif f[:2] == '-D':
+            tmp = string.split(f[2:], '=')
+            if len(tmp) == 1:
+                def_macros.append((tmp[0], None))
+            else:
+                def_macros.append(tuple(tmp))
+    return {'inc_dirs': inc_dirs, 'lib_dirs': lib_dirs, 'libs':libs,
+            'def_macros': def_macros, 'undef_macros': undef_macros}
+
+
+if __name__ == '__main__':
+    setup_compiler()
+
+    mpi_flags = get_mpi_flags()
+    mpi_flags['inc_dirs'].append(numpy.get_include())
+
+
+    # setting some extra compile flags for 64 bit architectures, utilizing
+    # distutils.sysconfig to check which compiler to use
+    if os.name == 'posix' and os.uname()[4] == 'x86_64':
+        #Extra flags for 64 bit architectures
+        extra_compile_args = ['-fPIC']
+    else:
+        extra_compile_args = None
+
+
+
+    setup(name='Pypar_extras',
+          version=__version__,
+          description='Extension to Pypar - Parallel Python',
+          long_description='Extension to Pypar - Parallel Python',
+          author=__author__,
+          author_email='stephen.roberts@anu.edu.au',
+          #url='http://sourceforge.net/projects/pypar',
+          package_dir = {'pypar_extras': ''}, # Use files in this dirctory
+          packages  = ['pypar_extras'],
+          ext_modules = [Extension('pypar_extras.mpiextras',
+                                   ['mpiextras.c'],
+                                   include_dirs=mpi_flags['inc_dirs'],
+                                   library_dirs=mpi_flags['lib_dirs'],
+                                   libraries=mpi_flags['libs'],
+                                   define_macros=mpi_flags['def_macros'],
+                                   undef_macros=mpi_flags['undef_macros'],
+                                   extra_compile_args=extra_compile_args)]
+         )
Index: /trunk/anuga_core/source/anuga_parallel/pypar_extras/test_pypar_ext.py
===================================================================
--- /trunk/anuga_core/source/anuga_parallel/pypar_extras/test_pypar_ext.py	(revision 8514)
+++ /trunk/anuga_core/source/anuga_parallel/pypar_extras/test_pypar_ext.py	(revision 8514)
@@ -0,0 +1,1257 @@
+# Test of MPI module 'pypar' for Python
+# 
+# Run as 
+#   python testpypar.py
+# or 
+#   mpirun -np 2 testpypar.py
+# (perhaps try number of processors more than 2)
+#
+# To verify bandwidth of your architecture please run pytiming (and ctiming) 
+#
+# OMN, GPC FEB 2002
+
+
+try:
+  import numpy
+except:
+  raise 'Module numpy must be present to run pypar'
+
+
+#print "Importing pypar"
+import pypar
+methods = dir(pypar)
+assert 'abort' in methods
+assert 'finalize' in methods
+assert 'get_processor_name' in methods
+assert 'time' in methods
+assert 'rank' in methods
+assert 'receive' in methods
+assert 'send' in methods
+assert 'broadcast' in methods
+assert 'size' in methods
+
+#print "Module pypar imported OK"
+#pypar.Barrier()
+
+
+myid = pypar.rank()
+numproc = pypar.size()
+node = pypar.get_processor_name()
+
+print 'I am processor %d of %d on node "%s"' %(myid, numproc, node)
+pypar.barrier()
+
+
+if numproc > 1:
+  # Test simple buffered communication (arrays, strings and general)
+  #
+  N = 17 #Number of elements
+  
+  if myid == 0:
+    # Integer arrays
+    #
+    A = numpy.array(range(N)).astype('i')
+    B = numpy.zeros(N).astype('i')
+
+    pypar.send(A, 1, use_buffer=True)
+    B = pypar.receive(numproc-1, buffer=B)
+
+    assert numpy.allclose(A, B)
+    print 'Buffered communication of numeric integer arrays OK'
+
+
+    # Long integer arrays
+    #
+    A = numpy.array(range(N)).astype('l')
+    B = numpy.zeros(N).astype('l')    
+
+    pypar.send(A, 1, use_buffer=True)
+    B = pypar.receive(numproc-1, buffer=B)
+    
+    assert numpy.allclose(A, B)
+    print 'Buffered communication of numeric long integer arrays OK'
+
+    # Real arrays
+    #
+    A = numpy.array(range(N)).astype('f')
+    B = numpy.zeros(N).astype('f')    
+
+    pypar.send(A, 1, use_buffer=True)
+    B = pypar.receive(numproc-1, buffer=B)
+    
+    assert numpy.allclose(A, B)    
+    print 'Buffered communication of numeric real arrays OK'
+
+    # Complex arrays
+    #
+    A = numpy.array(range(N)).astype('D')
+    A += 3j
+    B = numpy.zeros(N).astype('D')
+
+    pypar.send(A, 1, use_buffer=True)
+    B = pypar.receive(numproc-1, buffer=B)
+
+    assert numpy.allclose(A, B)    
+    print 'Buffered communication of numeric complex arrays OK'
+
+    # Strings (< 256 characters)
+    #
+    A = "and now to something completely different !"
+    B = " "*len(A)
+
+    pypar.send(A, 1, use_buffer=True)
+    B, status = pypar.receive(numproc-1, buffer=B, return_status=True)
+
+
+    assert A == B
+    print 'Buffered communication of strings OK'
+    
+    # A more general structure
+    #
+    A = ['ABC', (1,2,3.14), {8: 'Monty'}, numpy.array([13.45, 1.2])]
+    B = ['   ', (0,0,0.0) , {0: '     '}, numpy.zeros(2).astype('f')]    
+
+    pypar.send(A, 1, use_buffer=True)
+    B, status = pypar.receive(numproc-1, buffer=B, return_status=True)
+
+    OK = True
+    for i, a in enumerate(A):
+      b = B[i]
+
+      if type(a).__name__ == 'ndarray':
+        if not numpy.allclose(a, b):
+          OK = False
+          break
+      elif a != b:
+        OK = False
+        break
+
+    if OK is True:
+      print 'Buffered communication of general structures OK' 
+    else:
+      raise Exception    
+
+    
+  else:  
+    # Integers
+    #
+    X = numpy.zeros(N).astype('i')    
+
+    X = pypar.receive(myid-1, buffer=X)
+    pypar.send(X, (myid+1)%numproc, use_buffer=True)
+  
+    # Long integers
+    #
+    X = numpy.zeros(N).astype('l')
+    X = pypar.receive(myid-1, buffer=X)
+    pypar.send(X, (myid+1)%numproc, use_buffer=True)
+
+    # Floats
+    #
+    X = numpy.zeros(N).astype('f')
+    X = pypar.receive(myid-1, buffer=X)
+    pypar.send(X, (myid+1)%numproc, use_buffer=True)
+
+    # Complex
+    #
+    X = numpy.zeros(N).astype('D')
+    X = pypar.receive(myid-1, buffer=X)
+    pypar.send(X, (myid+1)%numproc, use_buffer=True)
+
+    # Strings
+    #
+    X = " "*256
+    X = pypar.receive(myid-1, buffer=X)
+    pypar.send(X.strip(), (myid+1)%numproc, use_buffer=True)
+
+    # General
+    #
+    X = ['   ', (0,0,0.0), {0: '     '}, numpy.zeros(2).astype('f')]
+    X = pypar.receive(myid-1, buffer=X)
+    pypar.send(X, (myid+1)%numproc, use_buffer=True)
+    
+
+
+
+
+  #Test (buffered communication of) multi dimensional arrays
+
+  M = 13  #Number of elements in dim 1
+  N = 17  #Number of elements in higher dims
+  
+  if myid == 0:
+    # 2D real arrays
+    #
+    A = numpy.array(range(M*N)).astype('f')
+    B = numpy.zeros(M*N).astype('f')
+
+    A = numpy.reshape(A, (M,N))
+    B = numpy.reshape(B, (M,N))    
+    
+    pypar.send(A, 1, use_buffer=True)
+    B = pypar.receive(numproc-1, buffer=B)
+    
+    assert numpy.allclose(A, B)    
+    print 'Buffered communication of 2D real arrays OK'
+
+    # 2D complex arrays
+    #
+    A = numpy.array(range(M*N)).astype('D')
+    B = numpy.zeros(M*N).astype('D')
+
+    A = numpy.reshape(A, (M,N))
+    B = numpy.reshape(B, (M,N))    
+
+    pypar.send(A, 1, use_buffer=True)
+    B = pypar.receive(numproc-1, buffer=B)
+    
+    assert numpy.allclose(A, B)    
+    print 'Buffered communication of 2D complex arrays OK'
+
+    # 3D real arrays
+    #
+    A = numpy.array(range(M*N*N)).astype('f')
+    B = numpy.zeros(M*N*N).astype('f')
+
+    A = numpy.reshape(A, (M,N,N))
+    B = numpy.reshape(B, (M,N,N))    
+    
+
+    pypar.send(A, 1, use_buffer=True)
+    B = pypar.receive(numproc-1, buffer=B)
+    
+    assert numpy.allclose(A, B)    
+    print 'Buffered communication of 3D real real arrays OK'
+
+    # 4D real arrays
+    #
+    A = numpy.array(range(M*N*N*M)).astype('f')
+    B = numpy.zeros(M*N*N*M).astype('f')
+
+    A = numpy.reshape(A, (M,N,N,M))
+    B = numpy.reshape(B, (M,N,N,M))    
+
+    pypar.send(A, 1, use_buffer=True)
+    B = pypar.receive(numproc-1, buffer=B)
+    
+    assert numpy.allclose(A, B)    
+    print 'Buffered communication of 4D real real arrays OK'
+
+    # 5D real arrays
+    #
+    A = numpy.array(range(M*N*2*N*M)).astype('f')
+    B = numpy.zeros(M*N*2*N*M).astype('f')
+
+    A = numpy.reshape(A, (M,N,2,N,M))
+    B = numpy.reshape(B, (M,N,2,N,M))    
+    
+    pypar.send(A, 1, use_buffer=True)
+    B = pypar.receive(numproc-1, buffer=B)
+    
+    assert numpy.allclose(A, B)    
+    print 'Buffered communication of 5D real real arrays OK'
+
+    # 5D double arrays
+    #
+    A = numpy.array(range(M*N*2*N*M)).astype('d')
+    B = numpy.zeros(M*N*2*N*M).astype('d')
+
+    A = numpy.reshape(A, (M,N,2,N,M))
+    B = numpy.reshape(B, (M,N,2,N,M))    
+    
+    pypar.send(A, 1, use_buffer=True)
+    B = pypar.receive(numproc-1, buffer=B)
+    
+    assert numpy.allclose(A, B)    
+    print 'Buffered communication of 5D double arrays OK'
+
+    # 5D complex arrays
+    #
+    A = numpy.array(range(M*N*2*N*M)).astype('D')
+    B = numpy.zeros(M*N*2*N*M).astype('D')
+
+    A = numpy.reshape(A, (M,N,2,N,M))
+    B = numpy.reshape(B, (M,N,2,N,M))    
+    
+    pypar.send(A, 1, use_buffer=True)
+    B = pypar.receive(numproc-1, buffer=B)
+    
+    assert numpy.allclose(A, B)    
+    print 'Buffered communication of 5D complex arrays OK'
+  else:  
+    # 2D real arrays
+    #
+    X = numpy.zeros(M*N).astype('f')
+    X = numpy.reshape(X, (M,N))
+    
+    X = pypar.receive(myid-1, buffer=X)
+    pypar.send(X, (myid+1)%numproc, use_buffer=True)
+  
+    # 2D complex arrays
+    #
+    X = numpy.zeros(M*N).astype('D')
+    X = numpy.reshape(X, (M,N))
+
+    X = pypar.receive(myid-1, buffer=X)
+    pypar.send(X, (myid+1)%numproc, use_buffer=True)
+  
+    # 3D real arrays
+    #
+    X = numpy.zeros(M*N*N).astype('f')
+    X = numpy.reshape(X, (M,N,N))
+    
+    X = pypar.receive(myid-1, buffer=X)
+    pypar.send(X, (myid+1)%numproc, use_buffer=True)
+
+    # 4D real arrays
+    #
+    X = numpy.zeros(M*N*N*M).astype('f')
+    X = numpy.reshape(X, (M,N,N,M))
+    
+    X = pypar.receive(myid-1, buffer=X)
+    pypar.send(X, (myid+1)%numproc, use_buffer=True)
+
+    # 5D real arrays
+    #
+    X = numpy.zeros(M*N*2*N*M).astype('f')
+    X = numpy.reshape(X, (M,N,2,N,M))
+    
+    X = pypar.receive(myid-1, buffer=X)
+    pypar.send(X, (myid+1)%numproc, use_buffer=True)
+
+    # 5D double arrays
+    #
+    X = numpy.zeros(M*N*2*N*M).astype('d')
+    X = numpy.reshape(X, (M,N,2,N,M))
+    
+    X = pypar.receive(myid-1, buffer=X)
+    pypar.send(X, (myid+1)%numproc, use_buffer=True)
+
+    # 5D complex arrays
+    #
+    X = numpy.zeros(M*N*2*N*M).astype('D')
+    X = numpy.reshape(X, (M,N,2,N,M))
+    
+    X = pypar.receive(myid-1, buffer=X)
+    pypar.send(X, (myid+1)%numproc, use_buffer=True)
+
+  # Test easy communication  - without buffers (arrays, strings and general)
+  #
+  N = 17 #Number of elements
+  
+  if myid == 0:
+    # Integer arrays
+    #
+    A = numpy.array(range(N))
+
+    pypar.send(A,1)
+    B = pypar.receive(numproc-1)
+    
+
+    assert numpy.allclose(A, B)
+    print "Simplified communication of numeric integer arrays OK"
+
+    # Long integer arrays
+    #
+    A = numpy.array(range(N)).astype('l')
+    pypar.send(A,1)
+    B=pypar.receive(numproc-1)
+    
+    assert numpy.allclose(A, B)    
+    print "Simplified communication of long integer real arrays OK"
+
+    # Real arrays
+    #
+    A = numpy.array(range(N)).astype('f')
+    pypar.send(A,1)
+    B=pypar.receive(numproc-1)
+    
+    assert numpy.allclose(A, B)    
+    print "Simplified communication of numeric real arrays OK"
+
+    # Complex arrays
+    #
+    A = numpy.array(range(N)).astype('D')
+    A += 3j
+    pypar.send(A,1)
+    B=pypar.receive(numproc-1)
+
+    assert numpy.allclose(A, B)    
+    print "Simplified communication of numeric complex arrays OK"
+
+
+    # Strings
+    #
+    A = "and now to something completely different !"
+    pypar.send(A,1)
+    B=pypar.receive(numproc-1)
+    
+    assert A == B
+    print "Simplified communication of strings OK"
+    
+    # A more general structure
+    #
+    A = ['ABC', (1,2,3.14), {8: 'Monty'}, numpy.array([13.45, 1.2])]
+    pypar.send(A,1)
+    B = pypar.receive(numproc-1)
+
+    OK = True
+    for i, a in enumerate(A):
+      b = B[i]
+
+      if type(a).__name__ == 'ndarray':
+        if not numpy.allclose(a, b):
+          OK = False
+          break
+      elif a != b:
+        OK = False
+        break
+
+    if OK is True:
+      print 'Simplified communication of general structures OK' 
+    else:
+      raise Exception    
+    
+  else:  
+    # Integers
+    #
+    X=pypar.receive(myid-1)  
+    pypar.send(X, (myid+1)%numproc)
+  
+    # Long integers
+    #
+    X=pypar.receive(myid-1)  
+    pypar.send(X, (myid+1)%numproc)    
+
+    # Floats
+    #
+    X=pypar.receive(myid-1)  
+    pypar.send(X, (myid+1)%numproc)    
+
+    # Complex
+    #
+    X=pypar.receive(myid-1)  
+    pypar.send(X, (myid+1)%numproc)    
+
+    # Strings
+    #
+    X=pypar.receive(myid-1)  
+    pypar.send(X, (myid+1)%numproc)    
+
+    # General
+    #
+    X = pypar.receive(myid-1)  
+    pypar.send(X, (myid+1)%numproc)    
+
+
+
+
+    
+  #Test (easy communication of) multi dimensional arrays
+
+  M = 13  #Number of elements in dim 1
+  N = 17  #Number of elements in higher dims
+  
+  if myid == 0:
+    # 2D real arrays
+    #
+    A = numpy.array(range(M*N)).astype('f')
+
+    A = numpy.reshape(A, (M,N))
+    
+    pypar.send(A,1)
+    B = pypar.receive(numproc-1)
+    
+    assert numpy.allclose(A, B)    
+    print "Simplified communication of 2D real arrays OK"
+
+    # 3D real arrays
+    #
+    A = numpy.array(range(M*N*N)).astype('f')
+    A = numpy.reshape(A, (M,N,N))
+    
+    pypar.send(A,1)
+    B=pypar.receive(numproc-1)
+    
+    assert numpy.allclose(A, B)    
+    print "Simplified communication of 3D real arrays OK"
+
+    # 4D real arrays
+    #
+    A = numpy.array(range(M*N*N*M)).astype('f')
+    A = numpy.reshape(A, (M,N,N,M))
+    
+    pypar.send(A,1)
+    B=pypar.receive(numproc-1)
+    
+    assert numpy.allclose(A, B)    
+    print "Simplified communication of 4D real arrays OK"
+
+    # 5D real arrays
+    #
+    A = numpy.array(range(M*N*2*N*M)).astype('f')
+    A = numpy.reshape(A, (M,N,2,N,M))
+    
+    pypar.send(A,1)
+    B=pypar.receive(numproc-1)
+    
+    assert numpy.allclose(A, B)    
+    print "Simplified communication of 5D real arrays OK"
+
+    # 5D double arrays
+    #
+    A = numpy.array(range(M*N*2*N*M)).astype('d')
+    A = numpy.reshape(A, (M,N,2,N,M))
+    
+    pypar.send(A,1)
+    B=pypar.receive(numproc-1)
+    
+    assert numpy.allclose(A, B)    
+    print "Simplified communication of 5D double arrays OK"
+
+    # 5D complex arrays
+    #
+    A = numpy.array(range(M*N*2*N*M)).astype('D')
+    A = numpy.reshape(A, (M,N,2,N,M))
+    
+    pypar.send(A,1)
+    B=pypar.receive(numproc-1)
+    
+    assert numpy.allclose(A, B)    
+    print "Simplified communication of 5D complex real arrays OK"
+
+  else:  
+    # 2D real arrays
+    #
+    
+    X = pypar.receive(myid-1)  
+    pypar.send(X, (myid+1)%numproc)
+  
+    # 3D real arrays
+    #
+    X=pypar.receive(myid-1)  
+    pypar.send(X, (myid+1)%numproc)
+
+    # 4D real arrays
+    #
+    X = pypar.receive(myid-1)  
+    pypar.send(X, (myid+1)%numproc)
+
+    # 5D real arrays
+    #
+    X=pypar.receive(myid-1)  
+    pypar.send(X, (myid+1)%numproc)
+
+    # 5D double arrays
+    #
+    X=pypar.receive(myid-1)  
+    pypar.send(X, (myid+1)%numproc)
+
+    # 5D complex arrays
+    #
+    X=pypar.receive(myid-1)  
+    pypar.send(X, (myid+1)%numproc)
+
+  # Test broadcast  - with buffers (arrays, strings and general)
+  #
+      
+  testString = ('test' + str(myid)).ljust(10)  #Buffers must have the same length on all procs!
+  pypar.broadcast(testString, 0)
+  assert testString.strip() == 'test0'
+  
+  testString = ('test' + str(myid)).ljust(10)  #Buffers must have the same length on all procs!
+  pypar.broadcast(testString, numproc-1)
+  assert testString.strip() == 'test' + str(numproc-1)
+  
+  if myid == 0:
+    print "Broadcast communication of strings OK"
+
+  ####################################################  
+  N = 17 #Number of elements
+  testArray = myid * numpy.array(range(N))
+  pypar.broadcast(testArray, 1)
+  assert numpy.allclose(testArray, 1 * testArray)
+  
+  if myid == 0:    
+    print "Broadcast communication of numeric integer array OK"
+
+
+  testArray = myid * numpy.array(range(N)).astype('f')
+  pypar.broadcast(testArray, 1)
+  assert numpy.allclose(testArray, 1 * testArray)
+  if myid == 0:
+    print "Broadcast communication of numeric real array OK"
+
+
+  M = 13
+  testArray = myid * numpy.array(range(M*N)).astype('f')
+  testArray = numpy.reshape(testArray, (M,N))  
+  pypar.broadcast(testArray, 1)
+  assert numpy.allclose(testArray, 1 * testArray)
+  if myid == 0:
+    print "Broadcast communication of 2D numeric real array OK"
+
+  testArray = myid * numpy.array(range(M*2*N)).astype('f')
+  testArray = numpy.reshape(testArray, (M,2,N))  
+  pypar.broadcast(testArray, 1)
+  assert numpy.allclose(testArray, 1 * testArray)
+  if myid == 0:
+    print "Broadcast communication of 3D numeric real array OK"
+
+  testArray = myid * numpy.array(range(M*2*N)).astype('D')
+  testArray = numpy.reshape(testArray, (M,2,N))  
+  pypar.broadcast(testArray, 1)
+  assert numpy.allclose(testArray, 1 * testArray)
+  if myid == 0:
+    print "Broadcast communication of 3D numeric complex array OK"
+
+    
+  A_x = ['ABC', myid, (1,2,3), {8: 'Monty'}, numpy.array([13.45, 1.2])]
+  A_1 = ['ABC',    1, (1,2,3), {8: 'Monty'}, numpy.array([13.45, 1.2])]
+  B = pypar.broadcast(A_x, 1)
+
+  OK = True
+  for i, a in enumerate(A_1):
+    b = B[i]
+
+    if type(a).__name__ == 'ndarray':
+      if not numpy.allclose(a, b):
+        OK = False
+        break
+    elif a != b:
+      OK = False
+      break
+
+  if OK is False:
+    raise Exception    
+    
+  if myid == 0:
+    print "Broadcast communication of general structures OK"
+  
+
+
+
+  # Test scatter  - with/without buffers (arrays, strings)
+  #
+  N = 16 # Number of elements
+
+  NP = N/numproc
+  
+  testString = 'ABCDEFGHIJKLMNOP'  #Length = 16
+  X = ' '*NP
+
+  pypar.scatter(testString, 0, buffer=X) # With buffer
+  Y = pypar.scatter(testString, 0) # With buffer automatically created
+
+  assert X==Y, 'X=%s, Y=%s' %(X,Y)
+  assert Y == testString[myid*NP:(myid+1)*NP]
+  assert X == testString[myid*NP:(myid+1)*NP]
+
+  if myid == 0:    
+    print "Scatter communication of strings OK"
+    
+
+  #Scatter Arrays
+  testArray = numpy.array(range(N)).astype('i')
+  X = numpy.zeros(NP).astype('i')
+
+  pypar.scatter(testArray, 0, buffer=X) # With buffer
+  Y = pypar.scatter(testArray, 0) # With buffer automatically created
+
+  assert numpy.allclose(X, Y)  
+  assert numpy.allclose(X, testArray[myid*NP:(myid+1)*NP])
+  assert numpy.allclose(Y, testArray[myid*NP:(myid+1)*NP])   
+
+  if myid == 0:
+    print "Scatter communication of numeric integer array OK"
+
+
+  testArray = numpy.array(range(N)).astype('f')
+  X = numpy.zeros(NP).astype('f')
+
+  pypar.scatter(testArray, 0, buffer=X) # With buffer
+  Y = pypar.scatter(testArray, 0) # With buffer automatically created
+
+  assert numpy.allclose(X, Y)  
+  assert numpy.allclose(X, testArray[myid*NP:(myid+1)*NP])
+  assert numpy.allclose(Y, testArray[myid*NP:(myid+1)*NP])   
+
+  if myid == 0:
+    print "Scatter communication of numeric real arrays OK"
+  #else:
+  #  print X, testArray, Y
+  #  assert numpy.allclose(X, Y)
+
+  testArray = numpy.array(range(N)).astype('D')
+  X = numpy.zeros(NP).astype('D')
+
+  pypar.scatter(testArray, 0, buffer=X) # With buffer
+  Y = pypar.scatter(testArray, 0) # With buffer automatically created
+
+  assert numpy.allclose(X, Y)  
+  assert numpy.allclose(X, testArray[myid*NP:(myid+1)*NP])
+  assert numpy.allclose(Y, testArray[myid*NP:(myid+1)*NP])   
+
+  if myid == 0:
+    print "Scatter communication of numeric complex array OK"
+
+
+###################################################
+#FIXME: 2D scatter doesn't work yet    
+#   M = 16
+#   N = 13
+#   MP = M/numproc
+  
+#   testArray = numpy.array(range(M*N)).astype('D')
+#   testArray = numpy.reshape(testArray, (M,N))    
+#   X = numpy.zeros(MP*N).astype('D')
+#   X = numpy.reshape(X, (MP,N))
+  
+#   raw_scatter(testArray, X, 0)
+#   Y = pypar.scatter(testArray, 0)
+  
+#   assert numpy.allclose(X, Y)  
+#   assert numpy.allclose(X, testArray[myid*MP:(myid+1)*MP,:])
+#   assert numpy.allclose(Y, testArray[myid*MP:(myid+1)*MP,:])   
+
+#   if myid == 0:
+#     print "Scatter communication of 2D numeric complex array OK"
+
+
+
+  # Test gather  - with/without buffers (arrays, strings)
+  #
+  N = 17 # Number of elements
+      
+  testString = 'AB'
+  X = '_'*(len(testString)*numproc) # Blanks caused errors when numproc >= 6 
+
+  pypar.gather(testString, 0, buffer=X) # Using buffer  
+  Y = pypar.gather(testString, 0) # Create buffer automatically
+  
+  if myid == 0:
+    assert X == 'AB' * numproc
+    assert Y == 'AB' * numproc
+    print "Gather communication of strings OK"
+  
+
+  testArray = numpy.array(range(N)).astype('i')
+  X = numpy.zeros(N*numproc).astype('i')
+
+  pypar.gather(testArray, 0, buffer=X) # Using buffer  
+  Y = pypar.gather(testArray, 0) # Create buffer automatically
+
+  if myid == 0:
+    for i in range(numproc):       
+      assert numpy.allclose(testArray, X[(i * N): ((i+1)*N)])
+
+    assert numpy.allclose(X, Y)
+    print "Gather communication of numeric integer array OK"
+    
+    
+  testArray = numpy.array(range(N)).astype('f')
+  X = numpy.zeros(N*numproc).astype('f')
+
+  pypar.gather(testArray, 0, buffer=X) # Using buffer  
+  Y = pypar.gather(testArray, 0) # Create buffer automatically
+
+  if myid == 0:
+    for i in range(numproc):       
+      assert numpy.allclose(testArray, X[(i * N): ((i+1)*N)])
+    assert numpy.allclose(X, Y)      
+    print "Gather communication of numeric real array OK"
+    
+  
+  testArray = numpy.array(range(N)).astype('D')
+  X = numpy.zeros(N*numproc).astype('D')
+
+  pypar.gather(testArray, 0, buffer=X) # Using buffer  
+  Y = pypar.gather(testArray, 0) # Create buffer automatically
+
+  if myid == 0:
+    for i in range(numproc):       
+      assert numpy.allclose(testArray, X[(i * N): ((i+1)*N)])
+    assert numpy.allclose(X, Y)      
+    print "Gather communication of numeric complex arrays OK"
+
+  M = 13  
+  testArray = numpy.array(range(M*N)).astype('D')
+  testArray = numpy.reshape(testArray, (M,N))
+  X = numpy.zeros(M*N*numproc).astype('D')
+  X = numpy.reshape(X, (M*numproc,N))
+  
+  pypar.gather(testArray, 0, buffer=X) # Using buffer  
+  Y = pypar.gather(testArray, 0) # Create buffer automatically
+
+  if myid == 0:
+    for i in range(numproc):       
+      assert numpy.allclose(testArray, X[(i * M): ((i+1)*M), :])
+    assert numpy.allclose(X, Y)      
+    print "Gather communication of 2D numeric complex arrays OK"
+  
+  ########################################################
+  # Test reduce
+  #######################################################
+  N = 17 #Number of elements
+  
+  # Create one (different) array on each processor
+  #    
+  testArray = numpy.array(range(N), 'i') * (myid+1)
+  X = numpy.zeros(N, 'i') # Buffer for results
+
+  pypar.reduce(testArray, pypar.SUM, 0, buffer=X)
+
+  if myid == 0:
+    Y = numpy.zeros(N, 'i')
+    for i in range(numproc):
+      Y = Y + numpy.array(range(N), 'i')*(i+1)    
+    assert numpy.allclose(X, Y)
+    print 'Buffered reduce using pypar.SUM OK'
+        
+  pypar.reduce(testArray, pypar.MAX, 0, buffer=X)
+  if myid == 0:
+    Y = numpy.array(range(N))*numproc
+    assert numpy.allclose(X, Y)
+    print 'Buffered reduce using pypar.MAX OK'
+
+  pypar.reduce(testArray, pypar.MIN, 0, buffer=X)
+  if myid == 0:
+    Y = numpy.array(range(N))
+    assert numpy.allclose(X, Y)
+    print 'Buffered reduce using pypar.MIN OK'
+    
+  if numproc <= 20:
+    testArray_float = testArray.astype('f')  #Make room for big results
+    X_float = X.astype('f')
+
+    pypar.reduce(testArray_float, pypar.PROD, 0, buffer=X_float)
+    if myid == 0:
+      Y = numpy.ones(N).astype('f')    
+      for i in range(numproc):
+        Y = Y*numpy.array(range(N))*(i+1)    
+      #print X_float
+      #print Y  
+      assert numpy.allclose(X_float, Y)
+      print 'Buffered reduce using pypar.PROD OK'
+  else:
+    if myid == 0:
+      print "Skipping product-reduce - try again with numproc < 20"    
+
+  pypar.reduce(testArray, pypar.LAND, 0, buffer=X)
+  if myid == 0:  
+    Y = numpy.ones(N).astype('i')    
+    for i in range(numproc):
+      Y = numpy.logical_and(Y, numpy.array(range(N)).astype('i')*(i+1))  
+    assert numpy.allclose(X, Y)
+    print 'Buffered reduce using pypar.LAND OK'
+    
+  pypar.reduce(testArray, pypar.BAND, 0, buffer=X)
+  if myid == 0:
+    Y = numpy.ones(N).astype('i')*255  #Neutral element for &   
+    for i in range(numproc):
+      Y = numpy.bitwise_and(Y, numpy.array(range(N))*(i+1))
+    assert numpy.allclose(X, Y)
+    print 'Buffered reduce using pypar.BAND OK'    
+
+  pypar.reduce(testArray, pypar.LOR, 0, buffer=X)
+  if myid == 0:  
+    Y = numpy.zeros(N).astype('i')    
+    for i in range(numproc):
+      Y = numpy.logical_or(Y, numpy.array(range(N)).astype('i')*(i+1))  
+    assert numpy.allclose(X, Y)
+    print 'Buffered reduce using pypar.LOR OK'
+  
+  pypar.reduce(testArray, pypar.BOR, 0, buffer=X)
+  if myid == 0:
+    Y = numpy.zeros(N).astype('i')   #Neutral element for |   
+    for i in range(numproc):
+      Y = numpy.bitwise_or(Y, numpy.array(range(N)).astype('i')*(i+1))
+    assert numpy.allclose(X, Y)
+    print 'Buffered reduce using pypar.BOR OK'    
+
+  pypar.reduce(testArray, pypar.LXOR, 0, buffer=X)
+  if myid == 0:  
+    Y = numpy.zeros(N).astype('i')    
+    for i in range(numproc):
+      Y = numpy.logical_xor(Y, numpy.array(range(N)).astype('i')*(i+1))  
+    assert numpy.allclose(X, Y)
+    print 'Buffered reduce using pypar.LXOR OK'
+
+  pypar.reduce(testArray, pypar.BXOR, 0, buffer=X)
+  if myid == 0:
+    Y = numpy.zeros(N).astype('i')   #Neutral element for xor ?   
+    for i in range(numproc):
+      Y = numpy.bitwise_xor(Y, numpy.array(range(N)).astype('i')*(i+1))
+    assert numpy.allclose(X, Y)
+    print 'Buffered reduce using pypar.BXOR OK'    
+
+
+
+
+
+  # FIXME: NOT YET SUPPORTED: pypar.MAXLOC, pypar.MINLOC, pypar.REPLACE
+
+  ########################################################
+  # Test allreduce
+  #######################################################
+  N = 17 #Number of elements
+
+  # Create one (different) array on each processor
+  #
+  testArray = numpy.array(range(N), 'i') * (myid+1)
+  X = numpy.zeros(N, 'i') # Buffer for results
+
+  pypar.allreduce(testArray, pypar.SUM, buffer=X)
+
+  if myid == 0:
+    Y = numpy.zeros(N, 'i')
+    for i in range(numproc):
+      Y = Y + numpy.array(range(N), 'i')*(i+1)
+    assert numpy.allclose(X, Y)
+    print 'Buffered allreduce using pypar.SUM OK'
+
+  pypar.allreduce(testArray, pypar.MAX, buffer=X)
+  if myid == 0:
+    Y = numpy.array(range(N))*numproc
+    assert numpy.allclose(X, Y)
+    print 'Buffered allreduce using pypar.MAX OK'
+
+  pypar.allreduce(testArray, pypar.MIN, buffer=X)
+  if myid == 0:
+    Y = numpy.array(range(N))
+    assert numpy.allclose(X, Y)
+    print 'Buffered allreduce using pypar.MIN OK'
+
+  if numproc <= 20:
+    testArray_float = testArray.astype('f')  #Make room for big results
+    X_float = X.astype('f')
+
+    pypar.allreduce(testArray_float, pypar.PROD, buffer=X_float)
+    if myid == 0:
+      Y = numpy.ones(N).astype('f')
+      for i in range(numproc):
+        Y = Y*numpy.array(range(N))*(i+1)
+      #print X_float
+      #print Y
+      assert numpy.allclose(X_float, Y)
+      print 'Buffered allreduce using pypar.PROD OK'
+  else:
+    if myid == 0:
+      print "Skipping product-allreduce - try again with numproc < 20"
+
+  pypar.allreduce(testArray, pypar.LAND, buffer=X)
+  if myid == 0:
+    Y = numpy.ones(N).astype('i')
+    for i in range(numproc):
+      Y = numpy.logical_and(Y, numpy.array(range(N)).astype('i')*(i+1))
+    assert numpy.allclose(X, Y)
+    print 'Buffered allreduce using pypar.LAND OK'
+
+  pypar.allreduce(testArray, pypar.BAND, buffer=X)
+  if myid == 0:
+    Y = numpy.ones(N).astype('i')*255  #Neutral element for &
+    for i in range(numproc):
+      Y = numpy.bitwise_and(Y, numpy.array(range(N))*(i+1))
+    assert numpy.allclose(X, Y)
+    print 'Buffered allreduce using pypar.BAND OK'
+
+  pypar.allreduce(testArray, pypar.LOR, buffer=X)
+  if myid == 0:
+    Y = numpy.zeros(N).astype('i')
+    for i in range(numproc):
+      Y = numpy.logical_or(Y, numpy.array(range(N)).astype('i')*(i+1))
+    assert numpy.allclose(X, Y)
+    print 'Buffered allreduce using pypar.LOR OK'
+
+  pypar.allreduce(testArray, pypar.BOR, buffer=X)
+  if myid == 0:
+    Y = numpy.zeros(N).astype('i')   #Neutral element for |
+    for i in range(numproc):
+      Y = numpy.bitwise_or(Y, numpy.array(range(N)).astype('i')*(i+1))
+    assert numpy.allclose(X, Y)
+    print 'Buffered allreduce using pypar.BOR OK'
+
+  pypar.allreduce(testArray, pypar.LXOR, buffer=X)
+  if myid == 0:
+    Y = numpy.zeros(N).astype('i')
+    for i in range(numproc):
+      Y = numpy.logical_xor(Y, numpy.array(range(N)).astype('i')*(i+1))
+    assert numpy.allclose(X, Y)
+    print 'Buffered allreduce using pypar.LXOR OK'
+
+  pypar.allreduce(testArray, pypar.BXOR, buffer=X)
+  if myid == 0:
+    Y = numpy.zeros(N).astype('i')   #Neutral element for xor ?
+    for i in range(numproc):
+      Y = numpy.bitwise_xor(Y, numpy.array(range(N)).astype('i')*(i+1))
+    assert numpy.allclose(X, Y)
+    print 'Buffered allreduce using pypar.BXOR OK'
+
+
+
+
+
+
+
+  # Test status block (simple communication)
+  N = 17  
+  if myid == 0:
+    # Integer arrays
+    #
+    A = numpy.array(range(N)).astype('i')
+
+    pypar.send(A,1)
+    B, status = pypar.receive(numproc-1, return_status = True)
+
+    repr(status)  # Check that status can be printed
+
+    sz = A.itemsize
+    assert numpy.allclose(A, B)
+    assert len(B) == status.length, 'Reported length == %d should be %d'\
+           %(status.length, len(B))
+    assert status.size == sz, 'Reported size == %d should be %d'\
+           %(status.size, sz)
+    assert status.tag == pypar.default_tag, 'Reported tag == %d should be %d'\
+           %(status.tag, pypar.default_tag)
+    assert status.error == 0, 'Status.error = %i' % status.error
+    assert status.source == numproc-1, 'Reported source == %d should be %d'\
+           %(status.source, numproc-1)
+
+    print "Status object (numeric integer arrays) OK"
+           
+    # Real arrays
+    #
+    A = numpy.array(range(N)).astype('f')
+    pypar.send(A,1)
+    B, status = pypar.receive(numproc-1, return_status = True)    
+
+    sz = A.itemsize    
+    assert numpy.allclose(A, B)
+    assert len(B) == status.length, 'Reported length == %d should be %d'\
+           %(status.length, len(B))
+    assert status.size == sz, 'Reported size == %d should be %d'\
+           %(status.size, sz)
+    assert status.tag == pypar.default_tag, 'Reported tag == %d should be %d'\
+           %(status.tag, pypar.default_tag)
+    assert status.error == 0
+    assert status.source == numproc-1, 'Reported source == %d should be %d'\
+           %(status.source, numproc-1)
+    
+    print "Status object (numeric real arrays) OK"
+
+    # Strings
+    #
+    A = "and now to something completely different !"
+    pypar.send(A,1)
+    B, status = pypar.receive(numproc-1, return_status = True)        
+
+    sz = 1 #Characters are one byte long
+    assert A == B
+    assert len(B) == status.length, 'Reported length == %d should be %d'\
+           %(status.length, len(B))
+    assert status.size == sz, 'Reported size == %d should be %d'\
+           %(status.size, sz)
+    assert status.tag == pypar.default_tag, 'Reported tag == %d should be %d'\
+           %(status.tag, pypar.default_tag)
+    assert status.error == 0
+    assert status.source == numproc-1, 'Reported source == %d should be %d'\
+           %(status.source, numproc-1)
+    
+    print "Status object (strings) OK"
+    
+    # A more general structure
+    #
+    A = ['ABC', (1,2,3.14), {8: 'Monty'}, numpy.array([13.45, 1.2])]
+    pypar.send(A,1)
+    B, status = pypar.receive(numproc-1, return_status = True)            
+
+    OK = True
+    for i, a in enumerate(A):
+      b = B[i]
+
+      if type(a).__name__ == 'ndarray':
+        if not numpy.allclose(a, b):
+          OK = False
+          break
+      elif a != b:
+        OK = False
+        break
+
+    if OK is True:
+      print 'Status object (more general structures) OK' 
+    else:
+      raise Exception    
+    
+
+    #Length is the number of characters needed to encode the structure
+    #Can't think of a test.
+
+    sz = 1
+    assert status.size == sz, 'Reported size == %d should be %d'\
+           %(status.size, sz)
+    assert status.tag == pypar.default_tag, 'Reported tag == %d should be %d'\
+           %(status.tag, pypar.default_tag)
+    assert status.error == 0
+    assert status.source == numproc-1, 'Reported source == %d should be %d'\
+           %(status.source, numproc-1)
+
+    print "Status object (general structures) OK"
+    
+  else:  
+    # Integers
+    #
+    X=pypar.receive(myid-1)  
+    pypar.send(X, (myid+1)%numproc)
+  
+    # Floats
+    #
+    X=pypar.receive(myid-1)  
+    pypar.send(X, (myid+1)%numproc)    
+
+    # Strings
+    #
+    X=pypar.receive(myid-1)  
+    pypar.send(X, (myid+1)%numproc)    
+
+    # General
+    #
+    X = pypar.receive(myid-1)  
+    pypar.send(X, (myid+1)%numproc)    
+
+
+
+
+
+
+
+
+  # Test status block (buffered communication)
+  N = 17 # Number of elements
+  if myid == 0:
+    # Integer arrays
+    #
+    A = numpy.array(range(N)).astype('i')
+    B = numpy.zeros(N).astype('i')    
+
+    pypar.send(A, 1, use_buffer=True)
+    B, status = pypar.receive(numproc-1, buffer=B, return_status=True)
+    
+    assert numpy.allclose(A, B)
+
+    sz = A.itemsize
+    assert numpy.allclose(A, B)
+    assert len(B) == status.length, 'Reported length == %d should be %d'\
+           %(status.length, len(B))
+    assert status.size == sz, 'Reported size == %d should be %d'\
+           %(status.size, sz)
+    assert status.tag == pypar.default_tag, 'Reported tag == %d should be %d'\
+           %(status.tag, pypar.default_tag)
+    assert status.error == 0
+    assert status.source == numproc-1, 'Reported source == %d should be %d'\
+           %(status.source, numproc-1)
+
+    print "Status object (buffered numeric integer arrays) OK"
+
+
+    # Real arrays
+    #
+    A = numpy.array(range(N)).astype('f')
+    B = numpy.zeros(N).astype('f')    
+    pypar.send(A, 1, use_buffer=True)
+    B, status = pypar.receive(numproc-1, buffer=B, return_status=True)
+    
+    assert numpy.allclose(A, B)
+    sz = A.itemsize
+    assert numpy.allclose(A, B)
+    assert len(B) == status.length, 'Reported length == %d should be %d'\
+           %(status.length, len(B))
+    assert status.size == sz, 'Reported size == %d should be %d'\
+           %(status.size, sz)
+    assert status.tag == pypar.default_tag, 'Reported tag == %d should be %d'\
+           %(status.tag, pypar.default_tag)
+    assert status.error == 0
+    assert status.source == numproc-1, 'Reported source == %d should be %d'\
+           %(status.source, numproc-1)
+
+    print "Status object (buffered numeric real arrays) OK"
+
+    # Strings (< 256 characters)
+    #
+    A = "and now to something completely different !"
+    B = " "*len(A)
+
+    pypar.send(A, 1, use_buffer=True)
+    B, status = pypar.receive(numproc-1, buffer=B, return_status=True)
+
+    sz = 1 #Characters are one byte long
+    assert A == B
+    assert len(B) == status.length, 'Reported length == %d should be %d'\
+           %(status.length, len(B))
+    assert status.size == sz, 'Reported size == %d should be %d'\
+           %(status.size, sz)
+    assert status.tag == pypar.default_tag, 'Reported tag == %d should be %d'\
+           %(status.tag, pypar.default_tag)
+    assert status.error == 0
+    assert status.source == numproc-1, 'Reported source == %d should be %d'\
+           %(status.source, numproc-1)
+
+    print "Status object (buffered strings) OK"
+    
+
+    
+    # A more general structure
+    #
+    A = ['ABC', (1,2,3.14), {8: 'Monty'}, numpy.array([13.45, 1.2])]
+    B = ['   ', (0,0,0.0), {0: '     '}, numpy.zeros(2).astype('f')]    
+    pypar.send(A, 1, use_buffer=True)
+    B, status = pypar.receive(numproc-1, buffer=B, return_status=True)
+
+    #assert A == B
+    sz = 1
+    assert status.size == sz, 'Reported size == %d should be %d'\
+           %(status.size, sz)
+    assert status.tag == pypar.default_tag, 'Reported tag == %d should be %d'\
+           %(status.tag, pypar.default_tag)
+    assert status.error == 0
+    assert status.source == numproc-1, 'Reported source == %d should be %d'\
+           %(status.source, numproc-1)
+
+    
+    print "Status object (buffered general structures) OK"
+    
+   
+  else:  
+    # Integers
+    #
+    X = numpy.zeros(N).astype('i')    
+
+
+    X = pypar.receive(myid-1, buffer=X)
+    pypar.send(X, (myid+1)%numproc, use_buffer=True)
+  
+    # Floats
+    #
+    X = numpy.zeros(N).astype('f')
+
+    X = pypar.receive(myid-1, buffer=X)
+    pypar.send(X, (myid+1)%numproc, use_buffer=True)
+
+    
+    # Strings
+    #
+    X = " "*256
+    X = pypar.receive(myid-1, buffer=X)
+    pypar.send(X.strip(), (myid+1)%numproc, use_buffer=True)
+
+    # General
+    #
+    X = ['   ', (0,0,0.0), {0: '     '}, numpy.zeros(2).astype('f')]
+    X = pypar.receive(myid-1, buffer=X)
+    pypar.send(X, (myid+1)%numproc, use_buffer=True)
+
+
+pypar.finalize()
+
+
+
+
