diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 4d46632e99820cc85b44076780142a6303686a10..1290ed98787377eea7d54f12d05b94eff0aeaf2f 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -18,14 +18,15 @@ workflow:
- if: $CI_COMMIT_BRANCH == $CI_DEFAULT_BRANCH
stages:
+ - compatibility
- build
- run
- test
-building_stage:
- stage: build
+compatibility_stage:
+ stage: compatibility
tags: ["np-tmcode"]
- allow_failure: false
+ allow_failure: true
artifacts:
paths:
- build/cluster/*
@@ -33,7 +34,6 @@ building_stage:
- build/testing/*
- build/trapping/*
- build/libnptm/*
- - doc/build/*
exclude:
- ".git*"
- ".git/**/*"
@@ -52,18 +52,44 @@ building_stage:
- echo "Running make with gnu compilers version 12..."
- make clean && BUILDDIR=$PWD/../build_gnu12 CXX=g++-12 FC=gfortran-12 make -j
- echo "Running make with flang version 16 and clang version 13..."
- - make clean && BUILDDIR=$PWD/../build_clang13-flang16 CXX="clang++-13 -stdlib=libc++" FC=flang-new-16 FCFLAGS=-O3 LDFLAGS=-L/usr/lib/llvm-16/lib make -j
+ - make clean && BUILDDIR=$PWD/../build_clang13-flang16 CXX="clang++-13 -stdlib=libstdc++ -I/usr/include/c++/12 -I/usr/include/x86_64-linux-gnu/c++/12" FC=flang-new-16 FCFLAGS=-O3 LDFLAGS="-L/usr/lib/llvm-16/lib -L/usr/lib/gcc/x86_64-linux-gnu/12" make -j
- echo "Running make with flang version 16 and clang version 14..."
- - make clean && BUILDDIR=$PWD/../build_clang14-flang16 CXX="clang++-14 -stdlib=libc++" FC=flang-new-16 FCFLAGS=-O3 LDFLAGS=-L/usr/lib/llvm-16/lib make -j
+ - make clean && BUILDDIR=$PWD/../build_clang14-flang16 CXX="clang++-14 -stdlib=libstdc++ -I/usr/include/c++/12 -I/usr/include/x86_64-linux-gnu/c++/12" FC=flang-new-16 FCFLAGS=-O3 LDFLAGS="-L/usr/lib/llvm-16/lib -L/usr/lib/gcc/x86_64-linux-gnu/12" make -j
- echo "Running make with flang version 16 and clang version 15..."
- - make clean && BUILDDIR=$PWD/../build_clang15-flang16 CXX="clang++-15 -stdlib=libc++" FC=flang-new-16 FCFLAGS=-O3 LDFLAGS=-L/usr/lib/llvm-16/lib make -j
+ - make clean && BUILDDIR=$PWD/../build_clang15-flang16 CXX="clang++-15 -stdlib=libstdc++ -I/usr/include/c++/12 -I/usr/include/x86_64-linux-gnu/c++/12" FC=flang-new-16 FCFLAGS=-O3 LDFLAGS="-L/usr/lib/llvm-16/lib -L/usr/lib/gcc/x86_64-linux-gnu/12" make -j
- echo "Running make with flang version 16 and clang version 16..."
- - make clean && BUILDDIR=$PWD/../build_clang16-flang16 CXX="clang++-16 -stdlib=libc++" FC=flang-new-16 FCFLAGS=-O3 LDFLAGS=-L/usr/lib/llvm-16/lib make -j
+ - make clean && BUILDDIR=$PWD/../build_clang16-flang16 CXX="clang++-16 -stdlib=libstdc++ -I/usr/include/c++/12 -I/usr/include/x86_64-linux-gnu/c++/12" FC=flang-new-16 FCFLAGS=-O3 LDFLAGS="-L/usr/lib/llvm-16/lib -L/usr/lib/gcc/x86_64-linux-gnu/12" make -j
- echo "Running make with Intel ifort and Intel icpx..."
- make clean && PATH=/opt/intel/oneapi/compiler/latest/bin:$PATH BUILDDIR=$PWD/../build_ifort-icpx CXX=icpx FC=ifort FCFLAGS="-O3 -diag-disable=10448" make -j
- echo "Running make with Intel ifx and Intel icpx..."
- make clean && LD_LIBRARY_PATH=/opt/intel/oneapi/compiler/latest/lib PATH=/opt/intel/oneapi/compiler/latest/bin:$PATH BUILDDIR=$PWD/../build_ifx-icpx CXX=icpx FC=ifx FCFLAGS=-O3 make -j
- - echo "Finally running make with default compilers..."
+
+building_stage:
+ stage: build
+ tags: ["np-tmcode"]
+ allow_failure: false
+ artifacts:
+ paths:
+ - build/cluster/*
+ - build/sphere/*
+ - build/testing/*
+ - build/trapping/*
+ - build/libnptm/*
+ - doc/build/*
+ exclude:
+ - ".git*"
+ - ".git/**/*"
+ expire_in: 2 hours
+ script:
+ # bash commands to be executed
+ - pwd
+ - hostname
+ - echo $CI_COMMIT_SHA
+ - echo $CI_COMMIT_BRANCH
+ - echo "Getting system info ..."
+ - cat /etc/os-release
+ - cd src
+ - echo "Running make with default compilers..."
- make clean && make -j
- make docs -j && make -C ../doc/build/latex -j
@@ -135,6 +161,9 @@ testing_stage:
- python3 ../../src/scripts/pycompare.py --ffile=$FFILE --cfile=c_OCLU --html
- echo "Testing cluster with 24 spheres"
- ./np_cluster ../../test_data/cluster/DEDFB_24 ../../test_data/cluster/DCLU_24 .
+ - echo "Comparing output of CLUSTER with 24 spheres"
+ - export FFILE=../../test_data/cluster/OCLU_24
+ - python3 ../../src/scripts/pycompare.py --ffile=$FFILE --cfile=c_OCLU --html
- echo "Checking consistency among legacy and HDF5 configuration files"
- ../testing/test_TEDF ../../test_data/cluster/DEDFB_24 c_TEDF c_TEDF.hd5
- echo "Checking consistency among legacy and HDF5 TM files"
diff --git a/COPYING b/COPYING
new file mode 100644
index 0000000000000000000000000000000000000000..f288702d2fa16d3cdf0035b15a9fcbc552cd88e7
--- /dev/null
+++ b/COPYING
@@ -0,0 +1,674 @@
+ GNU GENERAL PUBLIC LICENSE
+ Version 3, 29 June 2007
+
+ Copyright (C) 2007 Free Software Foundation, Inc. <https://fsf.org/>
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+ Preamble
+
+ The GNU General Public License is a free, copyleft license for
+software and other kinds of works.
+
+ The licenses for most software and other practical works are designed
+to take away your freedom to share and change the works. By contrast,
+the GNU General Public License is intended to guarantee your freedom to
+share and change all versions of a program--to make sure it remains free
+software for all its users. We, the Free Software Foundation, use the
+GNU General Public License for most of our software; it applies also to
+any other work released this way by its authors. You can apply it to
+your programs, too.
+
+ When we speak of free software, we are referring to freedom, not
+price. Our General Public Licenses are designed to make sure that you
+have the freedom to distribute copies of free software (and charge for
+them if you wish), that you receive source code or can get it if you
+want it, that you can change the software or use pieces of it in new
+free programs, and that you know you can do these things.
+
+ To protect your rights, we need to prevent others from denying you
+these rights or asking you to surrender the rights. Therefore, you have
+certain responsibilities if you distribute copies of the software, or if
+you modify it: responsibilities to respect the freedom of others.
+
+ For example, if you distribute copies of such a program, whether
+gratis or for a fee, you must pass on to the recipients the same
+freedoms that you received. You must make sure that they, too, receive
+or can get the source code. And you must show them these terms so they
+know their rights.
+
+ Developers that use the GNU GPL protect your rights with two steps:
+(1) assert copyright on the software, and (2) offer you this License
+giving you legal permission to copy, distribute and/or modify it.
+
+ For the developers' and authors' protection, the GPL clearly explains
+that there is no warranty for this free software. For both users' and
+authors' sake, the GPL requires that modified versions be marked as
+changed, so that their problems will not be attributed erroneously to
+authors of previous versions.
+
+ Some devices are designed to deny users access to install or run
+modified versions of the software inside them, although the manufacturer
+can do so. This is fundamentally incompatible with the aim of
+protecting users' freedom to change the software. The systematic
+pattern of such abuse occurs in the area of products for individuals to
+use, which is precisely where it is most unacceptable. Therefore, we
+have designed this version of the GPL to prohibit the practice for those
+products. If such problems arise substantially in other domains, we
+stand ready to extend this provision to those domains in future versions
+of the GPL, as needed to protect the freedom of users.
+
+ Finally, every program is threatened constantly by software patents.
+States should not allow patents to restrict development and use of
+software on general-purpose computers, but in those that do, we wish to
+avoid the special danger that patents applied to a free program could
+make it effectively proprietary. To prevent this, the GPL assures that
+patents cannot be used to render the program non-free.
+
+ The precise terms and conditions for copying, distribution and
+modification follow.
+
+ TERMS AND CONDITIONS
+
+ 0. Definitions.
+
+ "This License" refers to version 3 of the GNU General Public License.
+
+ "Copyright" also means copyright-like laws that apply to other kinds of
+works, such as semiconductor masks.
+
+ "The Program" refers to any copyrightable work licensed under this
+License. Each licensee is addressed as "you". "Licensees" and
+"recipients" may be individuals or organizations.
+
+ To "modify" a work means to copy from or adapt all or part of the work
+in a fashion requiring copyright permission, other than the making of an
+exact copy. The resulting work is called a "modified version" of the
+earlier work or a work "based on" the earlier work.
+
+ A "covered work" means either the unmodified Program or a work based
+on the Program.
+
+ To "propagate" a work means to do anything with it that, without
+permission, would make you directly or secondarily liable for
+infringement under applicable copyright law, except executing it on a
+computer or modifying a private copy. Propagation includes copying,
+distribution (with or without modification), making available to the
+public, and in some countries other activities as well.
+
+ To "convey" a work means any kind of propagation that enables other
+parties to make or receive copies. Mere interaction with a user through
+a computer network, with no transfer of a copy, is not conveying.
+
+ An interactive user interface displays "Appropriate Legal Notices"
+to the extent that it includes a convenient and prominently visible
+feature that (1) displays an appropriate copyright notice, and (2)
+tells the user that there is no warranty for the work (except to the
+extent that warranties are provided), that licensees may convey the
+work under this License, and how to view a copy of this License. If
+the interface presents a list of user commands or options, such as a
+menu, a prominent item in the list meets this criterion.
+
+ 1. Source Code.
+
+ The "source code" for a work means the preferred form of the work
+for making modifications to it. "Object code" means any non-source
+form of a work.
+
+ A "Standard Interface" means an interface that either is an official
+standard defined by a recognized standards body, or, in the case of
+interfaces specified for a particular programming language, one that
+is widely used among developers working in that language.
+
+ The "System Libraries" of an executable work include anything, other
+than the work as a whole, that (a) is included in the normal form of
+packaging a Major Component, but which is not part of that Major
+Component, and (b) serves only to enable use of the work with that
+Major Component, or to implement a Standard Interface for which an
+implementation is available to the public in source code form. A
+"Major Component", in this context, means a major essential component
+(kernel, window system, and so on) of the specific operating system
+(if any) on which the executable work runs, or a compiler used to
+produce the work, or an object code interpreter used to run it.
+
+ The "Corresponding Source" for a work in object code form means all
+the source code needed to generate, install, and (for an executable
+work) run the object code and to modify the work, including scripts to
+control those activities. However, it does not include the work's
+System Libraries, or general-purpose tools or generally available free
+programs which are used unmodified in performing those activities but
+which are not part of the work. For example, Corresponding Source
+includes interface definition files associated with source files for
+the work, and the source code for shared libraries and dynamically
+linked subprograms that the work is specifically designed to require,
+such as by intimate data communication or control flow between those
+subprograms and other parts of the work.
+
+ The Corresponding Source need not include anything that users
+can regenerate automatically from other parts of the Corresponding
+Source.
+
+ The Corresponding Source for a work in source code form is that
+same work.
+
+ 2. Basic Permissions.
+
+ All rights granted under this License are granted for the term of
+copyright on the Program, and are irrevocable provided the stated
+conditions are met. This License explicitly affirms your unlimited
+permission to run the unmodified Program. The output from running a
+covered work is covered by this License only if the output, given its
+content, constitutes a covered work. This License acknowledges your
+rights of fair use or other equivalent, as provided by copyright law.
+
+ You may make, run and propagate covered works that you do not
+convey, without conditions so long as your license otherwise remains
+in force. You may convey covered works to others for the sole purpose
+of having them make modifications exclusively for you, or provide you
+with facilities for running those works, provided that you comply with
+the terms of this License in conveying all material for which you do
+not control copyright. Those thus making or running the covered works
+for you must do so exclusively on your behalf, under your direction
+and control, on terms that prohibit them from making any copies of
+your copyrighted material outside their relationship with you.
+
+ Conveying under any other circumstances is permitted solely under
+the conditions stated below. Sublicensing is not allowed; section 10
+makes it unnecessary.
+
+ 3. Protecting Users' Legal Rights From Anti-Circumvention Law.
+
+ No covered work shall be deemed part of an effective technological
+measure under any applicable law fulfilling obligations under article
+11 of the WIPO copyright treaty adopted on 20 December 1996, or
+similar laws prohibiting or restricting circumvention of such
+measures.
+
+ When you convey a covered work, you waive any legal power to forbid
+circumvention of technological measures to the extent such circumvention
+is effected by exercising rights under this License with respect to
+the covered work, and you disclaim any intention to limit operation or
+modification of the work as a means of enforcing, against the work's
+users, your or third parties' legal rights to forbid circumvention of
+technological measures.
+
+ 4. Conveying Verbatim Copies.
+
+ You may convey verbatim copies of the Program's source code as you
+receive it, in any medium, provided that you conspicuously and
+appropriately publish on each copy an appropriate copyright notice;
+keep intact all notices stating that this License and any
+non-permissive terms added in accord with section 7 apply to the code;
+keep intact all notices of the absence of any warranty; and give all
+recipients a copy of this License along with the Program.
+
+ You may charge any price or no price for each copy that you convey,
+and you may offer support or warranty protection for a fee.
+
+ 5. Conveying Modified Source Versions.
+
+ You may convey a work based on the Program, or the modifications to
+produce it from the Program, in the form of source code under the
+terms of section 4, provided that you also meet all of these conditions:
+
+ a) The work must carry prominent notices stating that you modified
+ it, and giving a relevant date.
+
+ b) The work must carry prominent notices stating that it is
+ released under this License and any conditions added under section
+ 7. This requirement modifies the requirement in section 4 to
+ "keep intact all notices".
+
+ c) You must license the entire work, as a whole, under this
+ License to anyone who comes into possession of a copy. This
+ License will therefore apply, along with any applicable section 7
+ additional terms, to the whole of the work, and all its parts,
+ regardless of how they are packaged. This License gives no
+ permission to license the work in any other way, but it does not
+ invalidate such permission if you have separately received it.
+
+ d) If the work has interactive user interfaces, each must display
+ Appropriate Legal Notices; however, if the Program has interactive
+ interfaces that do not display Appropriate Legal Notices, your
+ work need not make them do so.
+
+ A compilation of a covered work with other separate and independent
+works, which are not by their nature extensions of the covered work,
+and which are not combined with it such as to form a larger program,
+in or on a volume of a storage or distribution medium, is called an
+"aggregate" if the compilation and its resulting copyright are not
+used to limit the access or legal rights of the compilation's users
+beyond what the individual works permit. Inclusion of a covered work
+in an aggregate does not cause this License to apply to the other
+parts of the aggregate.
+
+ 6. Conveying Non-Source Forms.
+
+ You may convey a covered work in object code form under the terms
+of sections 4 and 5, provided that you also convey the
+machine-readable Corresponding Source under the terms of this License,
+in one of these ways:
+
+ a) Convey the object code in, or embodied in, a physical product
+ (including a physical distribution medium), accompanied by the
+ Corresponding Source fixed on a durable physical medium
+ customarily used for software interchange.
+
+ b) Convey the object code in, or embodied in, a physical product
+ (including a physical distribution medium), accompanied by a
+ written offer, valid for at least three years and valid for as
+ long as you offer spare parts or customer support for that product
+ model, to give anyone who possesses the object code either (1) a
+ copy of the Corresponding Source for all the software in the
+ product that is covered by this License, on a durable physical
+ medium customarily used for software interchange, for a price no
+ more than your reasonable cost of physically performing this
+ conveying of source, or (2) access to copy the
+ Corresponding Source from a network server at no charge.
+
+ c) Convey individual copies of the object code with a copy of the
+ written offer to provide the Corresponding Source. This
+ alternative is allowed only occasionally and noncommercially, and
+ only if you received the object code with such an offer, in accord
+ with subsection 6b.
+
+ d) Convey the object code by offering access from a designated
+ place (gratis or for a charge), and offer equivalent access to the
+ Corresponding Source in the same way through the same place at no
+ further charge. You need not require recipients to copy the
+ Corresponding Source along with the object code. If the place to
+ copy the object code is a network server, the Corresponding Source
+ may be on a different server (operated by you or a third party)
+ that supports equivalent copying facilities, provided you maintain
+ clear directions next to the object code saying where to find the
+ Corresponding Source. Regardless of what server hosts the
+ Corresponding Source, you remain obligated to ensure that it is
+ available for as long as needed to satisfy these requirements.
+
+ e) Convey the object code using peer-to-peer transmission, provided
+ you inform other peers where the object code and Corresponding
+ Source of the work are being offered to the general public at no
+ charge under subsection 6d.
+
+ A separable portion of the object code, whose source code is excluded
+from the Corresponding Source as a System Library, need not be
+included in conveying the object code work.
+
+ A "User Product" is either (1) a "consumer product", which means any
+tangible personal property which is normally used for personal, family,
+or household purposes, or (2) anything designed or sold for incorporation
+into a dwelling. In determining whether a product is a consumer product,
+doubtful cases shall be resolved in favor of coverage. For a particular
+product received by a particular user, "normally used" refers to a
+typical or common use of that class of product, regardless of the status
+of the particular user or of the way in which the particular user
+actually uses, or expects or is expected to use, the product. A product
+is a consumer product regardless of whether the product has substantial
+commercial, industrial or non-consumer uses, unless such uses represent
+the only significant mode of use of the product.
+
+ "Installation Information" for a User Product means any methods,
+procedures, authorization keys, or other information required to install
+and execute modified versions of a covered work in that User Product from
+a modified version of its Corresponding Source. The information must
+suffice to ensure that the continued functioning of the modified object
+code is in no case prevented or interfered with solely because
+modification has been made.
+
+ If you convey an object code work under this section in, or with, or
+specifically for use in, a User Product, and the conveying occurs as
+part of a transaction in which the right of possession and use of the
+User Product is transferred to the recipient in perpetuity or for a
+fixed term (regardless of how the transaction is characterized), the
+Corresponding Source conveyed under this section must be accompanied
+by the Installation Information. But this requirement does not apply
+if neither you nor any third party retains the ability to install
+modified object code on the User Product (for example, the work has
+been installed in ROM).
+
+ The requirement to provide Installation Information does not include a
+requirement to continue to provide support service, warranty, or updates
+for a work that has been modified or installed by the recipient, or for
+the User Product in which it has been modified or installed. Access to a
+network may be denied when the modification itself materially and
+adversely affects the operation of the network or violates the rules and
+protocols for communication across the network.
+
+ Corresponding Source conveyed, and Installation Information provided,
+in accord with this section must be in a format that is publicly
+documented (and with an implementation available to the public in
+source code form), and must require no special password or key for
+unpacking, reading or copying.
+
+ 7. Additional Terms.
+
+ "Additional permissions" are terms that supplement the terms of this
+License by making exceptions from one or more of its conditions.
+Additional permissions that are applicable to the entire Program shall
+be treated as though they were included in this License, to the extent
+that they are valid under applicable law. If additional permissions
+apply only to part of the Program, that part may be used separately
+under those permissions, but the entire Program remains governed by
+this License without regard to the additional permissions.
+
+ When you convey a copy of a covered work, you may at your option
+remove any additional permissions from that copy, or from any part of
+it. (Additional permissions may be written to require their own
+removal in certain cases when you modify the work.) You may place
+additional permissions on material, added by you to a covered work,
+for which you have or can give appropriate copyright permission.
+
+ Notwithstanding any other provision of this License, for material you
+add to a covered work, you may (if authorized by the copyright holders of
+that material) supplement the terms of this License with terms:
+
+ a) Disclaiming warranty or limiting liability differently from the
+ terms of sections 15 and 16 of this License; or
+
+ b) Requiring preservation of specified reasonable legal notices or
+ author attributions in that material or in the Appropriate Legal
+ Notices displayed by works containing it; or
+
+ c) Prohibiting misrepresentation of the origin of that material, or
+ requiring that modified versions of such material be marked in
+ reasonable ways as different from the original version; or
+
+ d) Limiting the use for publicity purposes of names of licensors or
+ authors of the material; or
+
+ e) Declining to grant rights under trademark law for use of some
+ trade names, trademarks, or service marks; or
+
+ f) Requiring indemnification of licensors and authors of that
+ material by anyone who conveys the material (or modified versions of
+ it) with contractual assumptions of liability to the recipient, for
+ any liability that these contractual assumptions directly impose on
+ those licensors and authors.
+
+ All other non-permissive additional terms are considered "further
+restrictions" within the meaning of section 10. If the Program as you
+received it, or any part of it, contains a notice stating that it is
+governed by this License along with a term that is a further
+restriction, you may remove that term. If a license document contains
+a further restriction but permits relicensing or conveying under this
+License, you may add to a covered work material governed by the terms
+of that license document, provided that the further restriction does
+not survive such relicensing or conveying.
+
+ If you add terms to a covered work in accord with this section, you
+must place, in the relevant source files, a statement of the
+additional terms that apply to those files, or a notice indicating
+where to find the applicable terms.
+
+ Additional terms, permissive or non-permissive, may be stated in the
+form of a separately written license, or stated as exceptions;
+the above requirements apply either way.
+
+ 8. Termination.
+
+ You may not propagate or modify a covered work except as expressly
+provided under this License. Any attempt otherwise to propagate or
+modify it is void, and will automatically terminate your rights under
+this License (including any patent licenses granted under the third
+paragraph of section 11).
+
+ However, if you cease all violation of this License, then your
+license from a particular copyright holder is reinstated (a)
+provisionally, unless and until the copyright holder explicitly and
+finally terminates your license, and (b) permanently, if the copyright
+holder fails to notify you of the violation by some reasonable means
+prior to 60 days after the cessation.
+
+ Moreover, your license from a particular copyright holder is
+reinstated permanently if the copyright holder notifies you of the
+violation by some reasonable means, this is the first time you have
+received notice of violation of this License (for any work) from that
+copyright holder, and you cure the violation prior to 30 days after
+your receipt of the notice.
+
+ Termination of your rights under this section does not terminate the
+licenses of parties who have received copies or rights from you under
+this License. If your rights have been terminated and not permanently
+reinstated, you do not qualify to receive new licenses for the same
+material under section 10.
+
+ 9. Acceptance Not Required for Having Copies.
+
+ You are not required to accept this License in order to receive or
+run a copy of the Program. Ancillary propagation of a covered work
+occurring solely as a consequence of using peer-to-peer transmission
+to receive a copy likewise does not require acceptance. However,
+nothing other than this License grants you permission to propagate or
+modify any covered work. These actions infringe copyright if you do
+not accept this License. Therefore, by modifying or propagating a
+covered work, you indicate your acceptance of this License to do so.
+
+ 10. Automatic Licensing of Downstream Recipients.
+
+ Each time you convey a covered work, the recipient automatically
+receives a license from the original licensors, to run, modify and
+propagate that work, subject to this License. You are not responsible
+for enforcing compliance by third parties with this License.
+
+ An "entity transaction" is a transaction transferring control of an
+organization, or substantially all assets of one, or subdividing an
+organization, or merging organizations. If propagation of a covered
+work results from an entity transaction, each party to that
+transaction who receives a copy of the work also receives whatever
+licenses to the work the party's predecessor in interest had or could
+give under the previous paragraph, plus a right to possession of the
+Corresponding Source of the work from the predecessor in interest, if
+the predecessor has it or can get it with reasonable efforts.
+
+ You may not impose any further restrictions on the exercise of the
+rights granted or affirmed under this License. For example, you may
+not impose a license fee, royalty, or other charge for exercise of
+rights granted under this License, and you may not initiate litigation
+(including a cross-claim or counterclaim in a lawsuit) alleging that
+any patent claim is infringed by making, using, selling, offering for
+sale, or importing the Program or any portion of it.
+
+ 11. Patents.
+
+ A "contributor" is a copyright holder who authorizes use under this
+License of the Program or a work on which the Program is based. The
+work thus licensed is called the contributor's "contributor version".
+
+ A contributor's "essential patent claims" are all patent claims
+owned or controlled by the contributor, whether already acquired or
+hereafter acquired, that would be infringed by some manner, permitted
+by this License, of making, using, or selling its contributor version,
+but do not include claims that would be infringed only as a
+consequence of further modification of the contributor version. For
+purposes of this definition, "control" includes the right to grant
+patent sublicenses in a manner consistent with the requirements of
+this License.
+
+ Each contributor grants you a non-exclusive, worldwide, royalty-free
+patent license under the contributor's essential patent claims, to
+make, use, sell, offer for sale, import and otherwise run, modify and
+propagate the contents of its contributor version.
+
+ In the following three paragraphs, a "patent license" is any express
+agreement or commitment, however denominated, not to enforce a patent
+(such as an express permission to practice a patent or covenant not to
+sue for patent infringement). To "grant" such a patent license to a
+party means to make such an agreement or commitment not to enforce a
+patent against the party.
+
+ If you convey a covered work, knowingly relying on a patent license,
+and the Corresponding Source of the work is not available for anyone
+to copy, free of charge and under the terms of this License, through a
+publicly available network server or other readily accessible means,
+then you must either (1) cause the Corresponding Source to be so
+available, or (2) arrange to deprive yourself of the benefit of the
+patent license for this particular work, or (3) arrange, in a manner
+consistent with the requirements of this License, to extend the patent
+license to downstream recipients. "Knowingly relying" means you have
+actual knowledge that, but for the patent license, your conveying the
+covered work in a country, or your recipient's use of the covered work
+in a country, would infringe one or more identifiable patents in that
+country that you have reason to believe are valid.
+
+ If, pursuant to or in connection with a single transaction or
+arrangement, you convey, or propagate by procuring conveyance of, a
+covered work, and grant a patent license to some of the parties
+receiving the covered work authorizing them to use, propagate, modify
+or convey a specific copy of the covered work, then the patent license
+you grant is automatically extended to all recipients of the covered
+work and works based on it.
+
+ A patent license is "discriminatory" if it does not include within
+the scope of its coverage, prohibits the exercise of, or is
+conditioned on the non-exercise of one or more of the rights that are
+specifically granted under this License. You may not convey a covered
+work if you are a party to an arrangement with a third party that is
+in the business of distributing software, under which you make payment
+to the third party based on the extent of your activity of conveying
+the work, and under which the third party grants, to any of the
+parties who would receive the covered work from you, a discriminatory
+patent license (a) in connection with copies of the covered work
+conveyed by you (or copies made from those copies), or (b) primarily
+for and in connection with specific products or compilations that
+contain the covered work, unless you entered into that arrangement,
+or that patent license was granted, prior to 28 March 2007.
+
+ Nothing in this License shall be construed as excluding or limiting
+any implied license or other defenses to infringement that may
+otherwise be available to you under applicable patent law.
+
+ 12. No Surrender of Others' Freedom.
+
+ If conditions are imposed on you (whether by court order, agreement or
+otherwise) that contradict the conditions of this License, they do not
+excuse you from the conditions of this License. If you cannot convey a
+covered work so as to satisfy simultaneously your obligations under this
+License and any other pertinent obligations, then as a consequence you may
+not convey it at all. For example, if you agree to terms that obligate you
+to collect a royalty for further conveying from those to whom you convey
+the Program, the only way you could satisfy both those terms and this
+License would be to refrain entirely from conveying the Program.
+
+ 13. Use with the GNU Affero General Public License.
+
+ Notwithstanding any other provision of this License, you have
+permission to link or combine any covered work with a work licensed
+under version 3 of the GNU Affero General Public License into a single
+combined work, and to convey the resulting work. The terms of this
+License will continue to apply to the part which is the covered work,
+but the special requirements of the GNU Affero General Public License,
+section 13, concerning interaction through a network will apply to the
+combination as such.
+
+ 14. Revised Versions of this License.
+
+ The Free Software Foundation may publish revised and/or new versions of
+the GNU General Public License from time to time. Such new versions will
+be similar in spirit to the present version, but may differ in detail to
+address new problems or concerns.
+
+ Each version is given a distinguishing version number. If the
+Program specifies that a certain numbered version of the GNU General
+Public License "or any later version" applies to it, you have the
+option of following the terms and conditions either of that numbered
+version or of any later version published by the Free Software
+Foundation. If the Program does not specify a version number of the
+GNU General Public License, you may choose any version ever published
+by the Free Software Foundation.
+
+ If the Program specifies that a proxy can decide which future
+versions of the GNU General Public License can be used, that proxy's
+public statement of acceptance of a version permanently authorizes you
+to choose that version for the Program.
+
+ Later license versions may give you additional or different
+permissions. However, no additional obligations are imposed on any
+author or copyright holder as a result of your choosing to follow a
+later version.
+
+ 15. Disclaimer of Warranty.
+
+ THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY
+APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT
+HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY
+OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO,
+THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM
+IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF
+ALL NECESSARY SERVICING, REPAIR OR CORRECTION.
+
+ 16. Limitation of Liability.
+
+ IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS
+THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY
+GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE
+USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF
+DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD
+PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS),
+EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF
+SUCH DAMAGES.
+
+ 17. Interpretation of Sections 15 and 16.
+
+ If the disclaimer of warranty and limitation of liability provided
+above cannot be given local legal effect according to their terms,
+reviewing courts shall apply local law that most closely approximates
+an absolute waiver of all civil liability in connection with the
+Program, unless a warranty or assumption of liability accompanies a
+copy of the Program in return for a fee.
+
+ END OF TERMS AND CONDITIONS
+
+ 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 <https://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
+<https://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
+<https://www.gnu.org/licenses/why-not-lgpl.html>.
diff --git a/README.md b/README.md
index 15766ce3810b67aaa98dd64ac920233c06b0e1ee..cef20a5bc75eb48a494148946e1c6f438e2b86c6 100644
--- a/README.md
+++ b/README.md
@@ -7,3 +7,7 @@ The aim of the project, funded by PNRR-CNS, is to refactor the original, very ol
The current implementation offers a set of elementary tests to check that the original FORTRAN code can be compiled and executed on a limited set of pre-defined input data. The functionality of this initial stage can be verified by cloning the gitLab repository on a local machine and building the binaries from the `src` folder.
*NOTE:* The building process requires a working installation of the GNU Compiler Collection (`gcc`), of the GNU FORTRAN compiler (`gfortran`) and of the GNU `make` builder.
+
+# License
+
+This project is distributed under the terms of the *GNU General Public License V3.0 or later* (GPL-3.0-or-later). Check the contents of the *COPYING* file for further details.
\ No newline at end of file
diff --git a/build/README.md b/build/README.md
index 2596e698074726fbfbb035f1f95fe7a13d44335a..bd4fe07fd2a7db8b7b645a23a001180531071154 100644
--- a/build/README.md
+++ b/build/README.md
@@ -1,3 +1,5 @@
+Distributed under the terms of *GNU GPL-3.0-or-later*
+
# Folder instructions
This directory collects all the output of make builds.
diff --git a/containers/docker/Dockerfile b/containers/docker/Dockerfile
index 0d11a3db77d31a6b386bd70f968cd038efb82a6e..2f528d227e1881262a3297663a697b13d9ef5ffe 100644
--- a/containers/docker/Dockerfile
+++ b/containers/docker/Dockerfile
@@ -18,6 +18,8 @@ COPY --chown=root:root containers/docker/dockerstuff/intelcomps/oneapi-archive-k
COPY --chown=root:root containers/docker/dockerstuff/intelcomps/oneAPI.list /etc/apt/sources.list.d/
RUN apt update
RUN DEBIAN_FRONTEND=noninteractive apt -y install intel-oneapi-compiler-fortran intel-oneapi-compiler-dpcpp-cpp-and-cpp-classic intel-oneapi-compiler-dpcpp-cpp
+# install lapacke and its dependencies, both standard and the version with 64 bit integers
+RUN DEBIAN_FRONTEND=noninteractive apt -y install liblapacke-dev liblapacke64-dev libopenblas-dev libopenblas-openmp-dev libopenblas64-dev libopenblas64-openmp-dev
# install packages needed to run python scripts for checks
RUN DEBIAN_FRONTEND=noninteractive apt -y install python3 python-is-python3 python3-regex
# install packages needed to run doxygen to create html docs
@@ -51,7 +53,7 @@ FROM debian:bookworm-slim AS np-tmcode-run-minimal
WORKDIR /root
# install the strictly needed runtime libraries needed to run the executables
# and the python check scripts
-RUN DEBIAN_FRONTEND=noninteractive apt update && DEBIAN_FRONTEND=noninteractive apt upgrade && DEBIAN_FRONTEND=noninteractive apt -y install libgfortran5 libgcc-s1 libhdf5-103-1 libstdc++6 libssl3 libcurl4 libsz2 zlib1g libnghttp2-14 libidn2-0 librtmp1 libssh2-1 libpsl5 libgssapi-krb5-2 libldap-2.5-0 libzstd1 libbrotli1 libaec0 libunistring2 libgmp10 libkrb5-3 libk5crypto3 libcom-err2 libkrb5support0 libsasl2-2 libp11-kit0 libtasn1-6 libkeyutils1 libffi8 python3 python-is-python3 python3-regex hdf5-tools && rm -rf /var/lib/apt/lists/*
+RUN DEBIAN_FRONTEND=noninteractive apt update && DEBIAN_FRONTEND=noninteractive apt upgrade && DEBIAN_FRONTEND=noninteractive apt -y install libgfortran5 libgcc-s1 libhdf5-103-1 libstdc++6 libssl3 libcurl4 libsz2 zlib1g libnghttp2-14 libidn2-0 librtmp1 libssh2-1 libpsl5 libgssapi-krb5-2 libldap-2.5-0 libzstd1 libbrotli1 libaec0 libunistring2 libgmp10 libkrb5-3 libk5crypto3 libcom-err2 libkrb5support0 libsasl2-2 libp11-kit0 libtasn1-6 libkeyutils1 libffi8 liblapacke64 libopenblas64-0-openmp python3 python-is-python3 python3-regex hdf5-tools && rm -rf /var/lib/apt/lists/*
COPY --from=np-tmcode-run-dev /root /root
# remove everything which is not needed to run the codes
RUN cd /root/np-tmcode && find build -name "*.o" -exec rm -v \{\} \; && find build -name "*.gcno" -exec rm -v \{\} \; && cd src && rm -rvf cluster libnptm trapping include sphere Makefile make.inc README.md && cd .. && rm -rvf containers && cd doc && rm -rvf src && cd build/latex && rm -rvf *.tex *.out *.sty *.ind *.log *.toc *.ilg *.idx *.aux *.eps Makefile class*.pdf
diff --git a/containers/docker/README.md b/containers/docker/README.md
index c6bf1c9e9fc5bfbab6156a968259e6c0f66c386a..555a487d89a97d85acd42a412e535e3ba6e383ff 100644
--- a/containers/docker/README.md
+++ b/containers/docker/README.md
@@ -1,3 +1,5 @@
+Distributed under the terms of *GNU GPL-3.0-or-later*
+
# Docker support for the NP-TMcode project
This Dockerfile allows to create mainly:
diff --git a/containers/singularity/README.md b/containers/singularity/README.md
index 04cfdfe1108d63f93ebe707d6f35a5eb1956d580..20e9767c4df0e96126e41ef74e1b3b4b941dda42 100644
--- a/containers/singularity/README.md
+++ b/containers/singularity/README.md
@@ -1,3 +1,5 @@
+Distributed under the terms of *GNU GPL-3.0-or-later*
+
# Singularity support for the NP-TMcode project
This `np-tmcode-run.def` file allows to create a np-tmcode-run image, that contains only the pre-built executables, python test scripts, compiled documentation, and the minimal runtime to run them
diff --git a/containers/singularity/np-tmcode-run.def b/containers/singularity/np-tmcode-run.def
index da7025e0dd6899c8bf0698ec4615cb923bcc01ad..b9e3736ac27ea223c99de20775bbc7b739cc0a32 100644
--- a/containers/singularity/np-tmcode-run.def
+++ b/containers/singularity/np-tmcode-run.def
@@ -11,7 +11,7 @@ Stage: np-tmcode-run-dev
%post
apt update
apt -y upgrade
- apt -y install g++ gfortran make gcc-offload-nvptx libhdf5-dev
+ apt -y install g++ gfortran make gcc-offload-nvptx libhdf5-dev liblapacke-dev liblapacke64-dev libopenblas-dev libopenblas-openmp-dev libopenblas64-dev libopenblas64-openmp-dev
apt -y install python3 python-is-python3 python3-regex
apt -y install doxygen
apt -y install texlive-latex-base texlive-latex-recommended texlive-latex-extra texlive-font-utils
@@ -32,7 +32,7 @@ Stage: np-tmcode-run-minimal
%post
apt update
apt -y upgrade
- apt -y install libgfortran5 libgcc-s1 libhdf5-103-1 libstdc++6 libssl3 libcurl4 libsz2 zlib1g libnghttp2-14 libidn2-0 librtmp1 libssh2-1 libpsl5 libgssapi-krb5-2 libldap-2.5-0 libzstd1 libbrotli1 libaec0 libunistring2 libgmp10 libkrb5-3 libk5crypto3 libcom-err2 libkrb5support0 libsasl2-2 libp11-kit0 libtasn1-6 libkeyutils1 libffi8 python3 python-is-python3 python3-regex hdf5-tools
+ apt -y install libgfortran5 libgcc-s1 libhdf5-103-1 libstdc++6 libssl3 libcurl4 libsz2 zlib1g libnghttp2-14 libidn2-0 librtmp1 libssh2-1 libpsl5 libgssapi-krb5-2 libldap-2.5-0 libzstd1 libbrotli1 libaec0 libunistring2 libgmp10 libkrb5-3 libk5crypto3 libcom-err2 libkrb5support0 libsasl2-2 libp11-kit0 libtasn1-6 libkeyutils1 libffi8 liblapacke64 libopenblas64-0-openmp python3 python-is-python3 python3-regex hdf5-tools
rm -rf /var/lib/apt/lists/*
cd /usr/local/np-tmcode
find build -name "*.o" -exec rm -v \{\} \;
diff --git a/containers/singularity/np-tmcode-run.sif b/containers/singularity/np-tmcode-run.sif
index dd5129f356ed91f9477b7ee0c981b5b62cf6523d..d6fdf38eb04f927a4c9d721c1cab608f00b089c6 100755
Binary files a/containers/singularity/np-tmcode-run.sif and b/containers/singularity/np-tmcode-run.sif differ
diff --git a/doc/src/README.md b/doc/src/README.md
index 94056c1801324d3dcb85fffea980b27a26181a20..40ca9e988adb0865d58cb2011fb40b4bda0e3244 100644
--- a/doc/src/README.md
+++ b/doc/src/README.md
@@ -1,3 +1,5 @@
+Distributed under the terms of *GNU GPL-3.0-or-later*
+
# Folder instructions
This directory contains the material to build the project documentation with *doxygen*.
diff --git a/src/Makefile b/src/Makefile
index 61a13c6fe594c97d704dc0d1efaf62cc19a658cf..dd65769a321c70e4883df53b2e981f74936ee568 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -9,8 +9,8 @@ override BUILDDIR_NPTM=$(BUILDDIR)/libnptm
endif
ifndef LIBNPTM
# choose one of the two following lines, depending on whether a static or dynamic libnptm is wanted
-override LIBNPTM=$(BUILDDIR_NPTM)/libnptm.a
-#override LIBNPTM=$(BUILDDIR_NPTM)/libnptm.so
+#override LIBNPTM=$(BUILDDIR_NPTM)/libnptm.a
+override LIBNPTM=$(BUILDDIR_NPTM)/libnptm.so
endif
DOCSDIR=$(SRCDIR)/../doc
diff --git a/src/README.md b/src/README.md
index c345b3895299ffb296e3df970e74dcf92b3aca33..803470ee896cabebd92f61f297f971030136e1d4 100644
--- a/src/README.md
+++ b/src/README.md
@@ -1,3 +1,5 @@
+Distributed under the terms of *GNU GPL-3.0-or-later*
+
# Folder instructions
This directory collects the source code of the original programs and the development folders.
diff --git a/src/cluster/Makefile b/src/cluster/Makefile
index 22bc637a66fc67ed87f40b4355315a0e34ba87c8..c13c9680436ae9621e5b861ba2b7198138f20a1d 100644
--- a/src/cluster/Makefile
+++ b/src/cluster/Makefile
@@ -12,8 +12,8 @@ override BUILDDIR_NPTM=$(BUILDDIR)/libnptm
endif
ifndef LIBNPTM
# choose one of the two following lines, depending on whether a static or dynamic libnptm is wanted
-override LIBNPTM=$(BUILDDIR_NPTM)/libnptm.a
-#override LIBNPTM=$(BUILDDIR_NPTM)/libnpTm.so
+#override LIBNPTM=$(BUILDDIR_NPTM)/libnptm.a
+override LIBNPTM=$(BUILDDIR_NPTM)/libnptm.so
endif
include ../make.inc
@@ -46,7 +46,7 @@ $(BUILDDIR_CLU)/edfb_clu: $(OBJDIR) $(OBJDIR)/edfb_clu.o $(BUILDDIR_CLU)
# We put $(LIBNPTM) as an object to link in directly, so that it will be found at runtime even if it is a shared object library. May change in the future when we have an install: target
$(BUILDDIR_CLU)/np_cluster: $(OBJDIR) $(CXX_CLU_OBJS) $(BUILDDIR_CLU) $(LIBNPTM)
- $(CXX) $(CXXFLAGS) -o $(BUILDDIR_CLU)/np_cluster $(CXX_CLU_OBJS) $(LIBNPTM) $(CXXLDFLAGS)
+ $(CXX) $(CXXFLAGS) -o $(BUILDDIR_CLU)/np_cluster $(CXX_CLU_OBJS) $(LIBNPTM) $(CXXLDFLAGS)
clean:
rm -f $(F_CLU_OBJS) $(CXX_CLU_OBJS) $(CXX_CLU_DEBUG)
diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
index 4591b28aa6bda3f1483f0e0190af7bea94d934f2..e2debd9055ca0faf7a9711f360a7bd7d8f5637bd 100644
--- a/src/cluster/cluster.cpp
+++ b/src/cluster/cluster.cpp
@@ -1,13 +1,18 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file cluster.cpp
*
* \brief Implementation of the calculation for a cluster of spheres.
*/
-#include <complex>
#include <cstdio>
#include <exception>
#include <fstream>
#include <string>
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
+
#ifndef INCLUDE_ERRORS_H_
#include "../include/errors.h"
#endif
@@ -32,6 +37,21 @@
#include "../include/TransitionMatrix.h"
#endif
+#ifdef USE_LAPACK
+#ifdef USE_MKL
+#include <mkl_lapacke.h>
+#else
+#include <lapacke.h>
+#endif
+#ifndef INCLUDE_LAPACK_CALLS_H_
+#include "../include/lapack_calls.h"
+#endif
+#endif
+
+#ifndef INCLUDE_ALGEBRAIC_H_
+#include "../include/algebraic.h"
+#endif
+
using namespace std;
/*! \brief C++ implementation of CLU
@@ -66,7 +86,7 @@ void cluster(string config_file, string data_file, string output_path) {
if (sconf->number_of_spheres == gconf->number_of_spheres) {
// Shortcuts to variables stored in configuration objects
int nsph = gconf->number_of_spheres;
- int mxndm = gconf->mxndm;
+ lapack_int mxndm = (lapack_int)gconf->mxndm;
int inpol = gconf->in_pol;
int npnt = gconf->npnt;
int npntts = gconf->npntts;
@@ -117,49 +137,54 @@ void cluster(string config_file, string data_file, string output_path) {
C6 *c6 = new C6(c4->lmtpo);
FILE *output = fopen((output_path + "/c_OCLU").c_str(), "w");
int jer = 0, lcalc = 0;
- complex<double> arg(0.0, 0.0), ccsam(0.0, 0.0);
+ dcomplex arg = 0.0 + 0.0 * I;
+ dcomplex ccsam = 0.0 + 0.0 * I;
int configurations = 0;
for (int ci = 1; ci <= nsph; ci++) {
if (sconf->iog_vec[ci -1] >= ci) configurations++;
}
C2 *c2 = new C2(nsph, configurations, npnt, npntts);
- complex<double> **am = new complex<double>*[mxndm];
- for (int ai = 0; ai < mxndm; ai++) am[ai] = new complex<double>[mxndm]();
+ lapack_int ndit = 2 * nsph * c4->nlim;
+ dcomplex *am_vector = new dcomplex[ndit * ndit]();
+ dcomplex **am = new dcomplex*[ndit];
+ for (int ai = 0; ai < ndit; ai++) {
+ am[ai] = (am_vector + ai * ndit);
+ }
const int ndi = c4->nsph * c4->nlim;
C9 *c9 = new C9(ndi, c4->nlem, 2 * ndi, 2 * c4->nlem);
double *gaps = new double[nsph]();
double *tqev = new double[3]();
double *tqsv = new double[3]();
double **tqse, **tqss, **tqce, **tqcs;
- complex<double> **tqspe, **tqsps, **tqcpe, **tqcps;
+ dcomplex **tqspe, **tqsps, **tqcpe, **tqcps;
tqse = new double*[2];
- tqspe = new complex<double>*[2];
+ tqspe = new dcomplex*[2];
tqss = new double*[2];
- tqsps = new complex<double>*[2];
+ tqsps = new dcomplex*[2];
tqce = new double*[2];
- tqcpe = new complex<double>*[2];
+ tqcpe = new dcomplex*[2];
tqcs = new double*[2];
- tqcps = new complex<double>*[2];
+ tqcps = new dcomplex*[2];
for (int ti = 0; ti < 2; ti++) {
tqse[ti] = new double[nsph]();
- tqspe[ti] = new complex<double>[nsph]();
+ tqspe[ti] = new dcomplex[nsph]();
tqss[ti] = new double[nsph]();
- tqsps[ti] = new complex<double>[nsph]();
+ tqsps[ti] = new dcomplex[nsph]();
tqce[ti] = new double[3]();
- tqcpe[ti] = new complex<double>[3]();
+ tqcpe[ti] = new dcomplex[3]();
tqcs[ti] = new double[3]();
- tqcps[ti] = new complex<double>[3]();
+ tqcps[ti] = new dcomplex[3]();
}
double *gapv = new double[3]();
- complex<double> **gapp, **gappm;
+ dcomplex **gapp, **gappm;
double **gap, **gapm;
- gapp = new complex<double>*[3];
- gappm = new complex<double>*[3];
+ gapp = new dcomplex*[3];
+ gappm = new dcomplex*[3];
gap = new double*[3];
gapm = new double*[3];
for (int gi = 0; gi < 3; gi++) {
- gapp[gi] = new complex<double>[2]();
- gappm[gi] = new complex<double>[2]();
+ gapp[gi] = new dcomplex[2]();
+ gappm[gi] = new dcomplex[2]();
gap[gi] = new double[2]();
gapm[gi] = new double[2]();
}
@@ -192,7 +217,7 @@ void cluster(string config_file, string data_file, string output_path) {
double scan, cfmp, sfmp, cfsp, sfsp;
// End of global variables for CLU
fprintf(output, " READ(IR,*)NSPH,LI,LE,MXNDM,INPOL,NPNT,NPNTTS,IAVM,ISAM\n");
- fprintf(output, " %5d%5d%5d%5d%5d%5d%5d%5d%5d\n",
+ fprintf(output, " %5d%5d%5d%5ld%5d%5d%5d%5d%5d\n",
nsph, c4->li, c4->le, mxndm, inpol, npnt, npntts, iavm, isam
);
fprintf(output, " READ(IR,*)RXX(I),RYY(I),RZZ(I)\n");
@@ -264,6 +289,11 @@ void cluster(string config_file, string data_file, string output_path) {
string tppoan_name = output_path + "/c_TPPOAN";
tppoan.open(tppoan_name.c_str(), ios::out | ios::binary);
if (tppoan.is_open()) {
+#ifdef USE_LAPACK
+ printf("INFO: using LAPACK calls.\n");
+#else
+ printf("INFO: using fall-back lucin() calls.\n");
+#endif
tppoan.write(reinterpret_cast<char *>(&iavm), sizeof(int));
tppoan.write(reinterpret_cast<char *>(&isam), sizeof(int));
tppoan.write(reinterpret_cast<char *>(&inpol), sizeof(int));
@@ -331,8 +361,7 @@ void cluster(string config_file, string data_file, string output_path) {
if (jer != 0) break;
} // i132 loop
cms(am, c1, c1ao, c4, c6);
- int ndit = 2 * nsph * c4->nlim;
- lucin(am, mxndm, ndit, jer);
+ invert_matrix(am, ndit, jer, mxndm);
if (jer != 0) break; // jxi488 loop: goes to memory clean
ztm(am, c1, c1ao, c4, c6, c9);
if (sconf->idfc >= 0) {
@@ -354,9 +383,8 @@ void cluster(string config_file, string data_file, string output_path) {
// label 160
double cs0 = 0.25 * vk * vk * vk / acos(0.0);
double csch = 0.0, qschu = 0.0, pschu = 0.0, s0mag = 0.0;
- std::complex<double> s0(0.0, 0.0);
+ dcomplex s0 = 0.0 + 0.0 * I;
scr0(vk, exri, c1, c1ao, c3, c4);
- //printf("DEBUG: after SCR0 TFSAS = (%lE, %lE)\n", c3->tfsas.real(), c3->tfsas.imag());
double sqk = vk * vk * sconf->exdc;
aps(zpv, c4->li, nsph, c1, sqk, gaps);
rabas(inpol, c4->li, nsph, c1, tqse, tqspe, tqss, tqsps);
@@ -372,19 +400,19 @@ void cluster(string config_file, string data_file, string output_path) {
if (c1->nshl[i] != 1) {
fprintf(output, " SIZE=%15.7lE\n", c2->vsz[i]);
} else { // label 162
- fprintf(output, " SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE%15.7lE\n", c2->vsz[i], c2->vkt[i].real(), c2->vkt[i].imag());
+ fprintf(output, " SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE%15.7lE\n", c2->vsz[i], real(c2->vkt[i]), imag(c2->vkt[i]));
}
// label 164
fprintf(output, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n");
fprintf(output, " %14.7lE%15.7lE%15.7lE%15.7lE\n", c1->sscs[i], c1->sabs[i], c1->sexs[i], albeds);
fprintf(output, " ---- SCS/GS -- ABS/GS -- EXS/GS ---\n");
fprintf(output, " %14.7lE%15.7lE%15.7lE\n", c1->sqscs[i], c1->sqabs[i], c1->sqexs[i]);
- fprintf(output, " FSAS=%15.7lE%15.7lE\n", c1->fsas[i].real(), c1->fsas[i].imag());
+ fprintf(output, " FSAS=%15.7lE%15.7lE\n", real(c1->fsas[i]), imag(c1->fsas[i]));
csch = 2.0 * vk * sqsfi / c1->gcsv[i];
s0 = c1->fsas[i] * exri;
- qschu = s0.imag() * csch;
- pschu = s0.real() * csch;
- s0mag = abs(s0) * cs0;
+ qschu = imag(s0) * csch;
+ pschu = real(s0) * csch;
+ s0mag = cabs(s0) * cs0;
fprintf(output, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n", qschu, pschu, s0mag);
double rapr = c1->sexs[i] - gaps[i];
double cosav = gaps[i] / c1->sscs[i];
@@ -393,12 +421,12 @@ void cluster(string config_file, string data_file, string output_path) {
fprintf(output, " IPO= 2, TQEk=%15.7lE, TQSk=%15.7lE\n", tqse[1][i], tqss[1][i]);
}
} // i170 loop
- fprintf(output, " FSAT=%15.7lE%15.7lE\n", c3->tfsas.real(), c3->tfsas.imag());
+ fprintf(output, " FSAT=%15.7lE%15.7lE\n", real(c3->tfsas), imag(c3->tfsas));
csch = 2.0 * vk * sqsfi / c3->gcs;
s0 = c3->tfsas * exri;
- qschu = s0.imag() * csch;
- pschu = s0.real() * csch;
- s0mag = abs(s0) * cs0;
+ qschu = imag(s0) * csch;
+ pschu = real(s0) * csch;
+ s0mag = cabs(s0) * cs0;
fprintf(output, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n", qschu, pschu, s0mag);
tppoan.write(reinterpret_cast<char *>(&vk), sizeof(double));
pcrsm0(vk, exri, inpol, c1, c1ao, c4);
@@ -508,24 +536,24 @@ void cluster(string config_file, string data_file, string output_path) {
for (int i = 0; i < 2; i++) {
double value = c1ao->scscm[i];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = c1ao->scscpm[i].real();
+ value = real(c1ao->scscpm[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = c1ao->scscpm[i].imag();
+ value = imag(c1ao->scscpm[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
value = c1ao->ecscm[i];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = c1ao->ecscpm[i].real();
+ value = real(c1ao->ecscpm[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = c1ao->ecscpm[i].imag();
+ value = imag(c1ao->ecscpm[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
for (int i = 0; i < 3; i++) {
for (int j = 0; j < 2; j++) {
double value = gapm[i][j];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = gappm[i][j].real();
+ value = real(gappm[i][j]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = gappm[i][j].imag();
+ value = imag(gappm[i][j]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
}
@@ -546,12 +574,12 @@ void cluster(string config_file, string data_file, string output_path) {
double absrm = abssm / c3->acs;
double acsecs = c3->acs / c3->ecs;
if (acsecs >= -1.0e-6 && acsecs <= 1.0e-6) absrm = 1.0;
- complex<double> s0m = c1ao->fsacm[ilr210 - 1][ilr210 - 1] * exri;
- double qschum = s0m.imag() * csch;
- double pschum = s0m.real() * csch;
- double s0magm = abs(s0m) * cs0;
- double rfinrm = c1ao->fsacm[ilr210 - 1][ilr210 - 1].real() / c3->tfsas.real();
- double extcrm = c1ao->fsacm[ilr210 - 1][ilr210 - 1].imag() / c3->tfsas.imag();
+ dcomplex s0m = c1ao->fsacm[ilr210 - 1][ilr210 - 1] * exri;
+ double qschum = imag(s0m) * csch;
+ double pschum = real(s0m) * csch;
+ double s0magm = cabs(s0m) * cs0;
+ double rfinrm = real(c1ao->fsacm[ilr210 - 1][ilr210 - 1]) / real(c3->tfsas);
+ double extcrm = imag(c1ao->fsacm[ilr210 - 1][ilr210 - 1]) / imag(c3->tfsas);
if (inpol == 0) {
fprintf(output, " LIN %2d\n", ipol);
} else { // label 206
@@ -566,9 +594,9 @@ void cluster(string config_file, string data_file, string output_path) {
fprintf(output, " %14.7lE%15.7lE%15.7lE\n", scarm, absrm, extrm);
fprintf(
output, " FSAC(%1d,%1d)=%15.7lE%15.7lE FSAC(%1d,%1d)=%15.7lE%15.7lE\n",
- ilr210, ilr210, c1ao->fsacm[ilr210 - 1][ilr210 - 1].real(),
- c1ao->fsacm[ilr210 - 1][ilr210 - 1].imag(), jlr, ilr210,
- c1ao->fsacm[jlr - 1][ilr210 - 1].real(), c1ao->fsacm[jlr - 1][ilr210 - 1].imag()
+ ilr210, ilr210, real(c1ao->fsacm[ilr210 - 1][ilr210 - 1]),
+ imag(c1ao->fsacm[ilr210 - 1][ilr210 - 1]), jlr, ilr210,
+ real(c1ao->fsacm[jlr - 1][ilr210 - 1]), imag(c1ao->fsacm[jlr - 1][ilr210 - 1])
);
fprintf(
output, " RE(FSAC(%1d,%1d))/RE(TFSAS)=%15.7lE, IM(FSAC(%1d,%1d))/IM(TFSAS)=%15.7lE\n",
@@ -581,8 +609,8 @@ void cluster(string config_file, string data_file, string output_path) {
fprintf(output, " COSAV=%15.7lE, RAPRS=%15.7lE\n", cosav, rapr);
fprintf(output, " Fk=%15.7lE\n", fz);
} // ilr210 loop
- double rmbrif = (c1ao->fsacm[0][0].real() - c1ao->fsacm[1][1].real()) / c1ao->fsacm[0][0].real();
- double rmdchr = (c1ao->fsacm[0][0].imag() - c1ao->fsacm[1][1].imag()) / c1ao->fsacm[0][0].imag();
+ double rmbrif = (real(c1ao->fsacm[0][0]) - real(c1ao->fsacm[1][1])) / real(c1ao->fsacm[0][0]);
+ double rmdchr = (imag(c1ao->fsacm[0][0]) - imag(c1ao->fsacm[1][1])) / imag(c1ao->fsacm[0][0]);
fprintf(output, " (RE(FSAC(1,1))-RE(FSAC(2,2)))/RE(FSAC(1,1))=%15.7lE\n", rmbrif);
fprintf(output, " (IM(FSAC(1,1))-IM(FSAC(2,2)))/IM(FSAC(1,1))=%15.7lE\n", rmdchr);
}
@@ -612,13 +640,13 @@ void cluster(string config_file, string data_file, string output_path) {
fprintf(output, " SPHERE %2d\n", i226);
fprintf(
output, " SAS(1,1)=%15.7lE%15.7lE, SAS(2,1)=%15.7lE%15.7lE\n",
- c1->sas[i226 - 1][0][0].real(), c1->sas[i226 - 1][0][0].imag(),
- c1->sas[i226 - 1][1][0].real(), c1->sas[i226 - 1][1][0].imag()
+ real(c1->sas[i226 - 1][0][0]), imag(c1->sas[i226 - 1][0][0]),
+ real(c1->sas[i226 - 1][1][0]), imag(c1->sas[i226 - 1][1][0])
);
fprintf(
output, " SAS(1,2)=%15.7lE%15.7lE, SAS(2,2)=%15.7lE%15.7lE\n",
- c1->sas[i226 - 1][0][1].real(), c1->sas[i226 - 1][0][1].imag(),
- c1->sas[i226 - 1][1][1].real(), c1->sas[i226 - 1][1][1].imag()
+ real(c1->sas[i226 - 1][0][1]), imag(c1->sas[i226 - 1][0][1]),
+ real(c1->sas[i226 - 1][1][1]), imag(c1->sas[i226 - 1][1][1])
);
for (int j225 = 0; j225 < 16; j225++) { // QUESTION: check that 16 is a fixed dimension
c1ao->vint[j225] = c1ao->vints[i226 - 1][j225];
@@ -642,13 +670,13 @@ void cluster(string config_file, string data_file, string output_path) {
} // i226 loop
fprintf(
output, " SAT(1,1)=%15.7lE%15.7lE, SAT(2,1)=%15.7lE%15.7lE\n",
- c3->tsas[0][0].real(), c3->tsas[0][0].imag(),
- c3->tsas[1][0].real(), c3->tsas[1][0].imag()
+ real(c3->tsas[0][0]), imag(c3->tsas[0][0]),
+ real(c3->tsas[1][0]), imag(c3->tsas[1][0])
);
fprintf(
output, " SAT(1,2)=%15.7lE%15.7lE, SAT(2,2)=%15.7lE%15.7lE\n",
- c3->tsas[0][1].real(), c3->tsas[0][1].imag(),
- c3->tsas[1][1].real(), c3->tsas[1][1].imag()
+ real(c3->tsas[0][1]), imag(c3->tsas[0][1]),
+ real(c3->tsas[1][1]), imag(c3->tsas[1][1])
);
fprintf(output, " CLUSTER\n");
pcros(vk, exri, c1, c1ao, c4);
@@ -666,24 +694,24 @@ void cluster(string config_file, string data_file, string output_path) {
for (int i = 0; i < 2; i++) {
double value = c1ao->scsc[i];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = c1ao->scscp[i].real();
+ value = real(c1ao->scscp[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = c1ao->scscp[i].imag();
+ value = imag(c1ao->scscp[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
value = c1ao->ecsc[i];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = c1ao->ecscp[i].real();
+ value = real(c1ao->ecscp[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = c1ao->ecscp[i].imag();
+ value = imag(c1ao->ecscp[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
for (int i = 0; i < 3; i++) {
for (int j = 0; j < 2; j++) {
double value = gap[i][j];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = gapp[i][j].real();
+ value = real(gapp[i][j]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = gapp[i][j].imag();
+ value = imag(gapp[i][j]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
}
@@ -691,9 +719,9 @@ void cluster(string config_file, string data_file, string output_path) {
for (int j = 0; j < 3; j++) {
double value = tqce[i][j];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = tqcpe[i][j].real();
+ value = real(tqcpe[i][j]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = tqcpe[i][j].imag();
+ value = imag(tqcpe[i][j]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
}
@@ -701,9 +729,9 @@ void cluster(string config_file, string data_file, string output_path) {
for (int j = 0; j < 3; j++) {
double value = tqcs[i][j];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = tqcps[i][j].real();
+ value = real(tqcps[i][j]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = tqcps[i][j].imag();
+ value = imag(tqcps[i][j]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
}
@@ -718,9 +746,9 @@ void cluster(string config_file, string data_file, string output_path) {
}
// label 254
for (int i = 0; i < 16; i++) {
- double value = c1ao->vint[i].real();
+ double value = real(c1ao->vint[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = c1ao->vint[i].imag();
+ value = imag(c1ao->vint[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
for (int i = 0; i < 4; i++) {
@@ -746,11 +774,11 @@ void cluster(string config_file, string data_file, string output_path) {
double ratio = c3->acs / c3->ecs;
if (ratio < -1.0e-6 || ratio > 1.0e-6) absrat = abssec / c3->acs;
s0 = c1ao->fsac[ilr290 - 1][ilr290 - 1] * exri;
- double qschu = s0.imag() * csch;
- double pschu = s0.real() * csch;
- s0mag = abs(s0) * cs0;
- double refinr = c1ao->fsac[ilr290 - 1][ilr290 - 1].real() / c3->tfsas.real();
- double extcor = c1ao->fsac[ilr290 - 1][ilr290 - 1].imag() / c3->tfsas.imag();
+ double qschu = imag(s0) * csch;
+ double pschu = real(s0) * csch;
+ s0mag = cabs(s0) * cs0;
+ double refinr = real(c1ao->fsac[ilr290 - 1][ilr290 - 1]) / real(c3->tfsas);
+ double extcor = imag(c1ao->fsac[ilr290 - 1][ilr290 - 1]) / imag(c3->tfsas);
if (inpol == 0) {
fprintf(output, " LIN %2d\n", ipol);
} else { // label 273
@@ -774,13 +802,13 @@ void cluster(string config_file, string data_file, string output_path) {
);
fprintf(
output, " FSAC(%1d,%1d)=%15.7lE%15.7lE FSAC(%1d,%1d)=%15.7lE%15.7lE\n",
- ilr290, ilr290, c1ao->fsac[ilr290 - 1][ilr290 - 1].real(), c1ao->fsac[ilr290 - 1][ilr290 - 1].imag(),
- jlr, ilr290, c1ao->fsac[jlr - 1][ilr290 - 1].real(), c1ao->fsac[jlr - 1][ilr290 - 1].imag()
+ ilr290, ilr290, real(c1ao->fsac[ilr290 - 1][ilr290 - 1]), imag(c1ao->fsac[ilr290 - 1][ilr290 - 1]),
+ jlr, ilr290, real(c1ao->fsac[jlr - 1][ilr290 - 1]), imag(c1ao->fsac[jlr - 1][ilr290 - 1])
);
fprintf(
output, " SAC(%1d,%1d)=%15.7lE%15.7lE SAC(%1d,%1d)=%15.7lE%15.7lE\n",
- ilr290, ilr290, c1ao->sac[ilr290 - 1][ilr290 - 1].real(), c1ao->sac[ilr290 - 1][ilr290 - 1].imag(),
- jlr, ilr290, c1ao->sac[jlr - 1][ilr290 - 1].real(), c1ao->sac[jlr - 1][ilr290 - 1].imag()
+ ilr290, ilr290, real(c1ao->sac[ilr290 - 1][ilr290 - 1]), imag(c1ao->sac[ilr290 - 1][ilr290 - 1]),
+ jlr, ilr290, real(c1ao->sac[jlr - 1][ilr290 - 1]), imag(c1ao->sac[jlr - 1][ilr290 - 1])
);
fprintf(
output, " RE(FSAC(%1d,%1d))/RE(TFSAS)=%15.7lE, IM(FSAC(%1d,%1d))/IM(TFSAS)=%15.7lE\n",
@@ -822,8 +850,8 @@ void cluster(string config_file, string data_file, string output_path) {
);
}
} //ilr290 loop
- double rbirif = (c1ao->fsac[0][0].real() - c1ao->fsac[1][1].real()) / c1ao->fsac[0][0].real();
- double rdichr = (c1ao->fsac[0][0].imag() - c1ao->fsac[1][1].imag()) / c1ao->fsac[0][0].imag();
+ double rbirif = (real(c1ao->fsac[0][0]) - real(c1ao->fsac[1][1])) / real(c1ao->fsac[0][0]);
+ double rdichr = (imag(c1ao->fsac[0][0]) - imag(c1ao->fsac[1][1])) / imag(c1ao->fsac[0][0]);
fprintf(output, " (RE(FSAC(1,1))-RE(FSAC(2,2)))/RE(FSAC(1,1))=%15.7lE\n", rbirif);
fprintf(output, " (IM(FSAC(1,1))-IM(FSAC(2,2)))/IM(FSAC(1,1))=%15.7lE\n", rdichr);
fprintf(output, " MULC\n");
@@ -845,9 +873,9 @@ void cluster(string config_file, string data_file, string output_path) {
// Some implicit loops writing to binary.
for (int i = 0; i < 16; i++) {
double value;
- value = c1ao->vintm[i].real();
+ value = real(c1ao->vintm[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = c1ao->vintm[i].imag();
+ value = imag(c1ao->vintm[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
for (int i = 0; i < 4; i++) {
@@ -910,8 +938,9 @@ void cluster(string config_file, string data_file, string output_path) {
delete[] zpv[zi];
}
delete[] zpv;
- for (int ai = mxndm - 1; ai > -1; ai--) delete[] am[ai];
+ delete[] am_vector;
delete[] am;
+ //delete[] tam;
delete[] gaps;
for (int ti = 1; ti > -1; ti--) {
delete[] tqse[ti];
diff --git a/src/cluster/np_cluster.cpp b/src/cluster/np_cluster.cpp
index 73caa518f9d330b8ad4114a30f123ba320035756..66fb783bc2e9d51e32b5dbc6fb0c033793348936 100644
--- a/src/cluster/np_cluster.cpp
+++ b/src/cluster/np_cluster.cpp
@@ -1,3 +1,5 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file np_cluster.cpp
*
* \brief Cluster of spheres scattering problem handler.
@@ -17,9 +19,12 @@
*/
#include <cstdio>
-#include <complex>
#include <string>
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
+
#ifndef INCLUDE_CONFIGURATION_H_
#include "../include/Configuration.h"
#endif
diff --git a/src/include/Commons.h b/src/include/Commons.h
index 6ad2a0110e6eb4ebae592402e9df43e9cb182687..e954544404d3ecc19c88377897b18c2d7e0ea1d8 100644
--- a/src/include/Commons.h
+++ b/src/include/Commons.h
@@ -1,3 +1,5 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file Commons.h
*
* \brief C++ porting of common data structures.
@@ -42,15 +44,15 @@ protected:
int nlmmt;
public:
//! \brief QUESTION: definition?
- std::complex<double> **rmi;
+ dcomplex **rmi;
//! \brief QUESTION: definition?
- std::complex<double> **rei;
+ dcomplex **rei;
//! \brief QUESTION: definition?
- std::complex<double> **w;
+ dcomplex **w;
//! \brief QUESTION: definition?
- std::complex<double> *fsas;
+ dcomplex *fsas;
//! \brief QUESTION: definition?
- std::complex<double> **vints;
+ dcomplex **vints;
//! \brief QUESTION: definition?
double *sscs;
//! \brief QUESTION: definition?
@@ -80,7 +82,7 @@ public:
//! \brief Vector of numbers of spherical layers.
int *nshl;
//! \brief QUESTION: definition?
- std::complex<double> ***sas;
+ dcomplex ***sas;
/*! \brief C1 instance constructor.
*
@@ -105,13 +107,13 @@ class C2 {
int nhspo;
public:
//! \brief QUESTION: definition?
- std::complex<double> *ris;
+ dcomplex *ris;
//! \brief QUESTION: definition?
- std::complex<double> *dlri;
+ dcomplex *dlri;
//! \brief QUESTION: definition?
- std::complex<double> *dc0;
+ dcomplex *dc0;
//! \brief QUESTION: definition?
- std::complex<double> *vkt;
+ dcomplex *vkt;
//! Vector of scaled sizes. QUESTION: correct?
double *vsz;
@@ -133,9 +135,9 @@ public:
class C3 {
public:
//! \brief QUESTION: definition?
- std::complex<double> tfsas;
+ dcomplex tfsas;
//! \brief QUESTION: definition?
- std::complex<double> **tsas;
+ dcomplex **tsas;
//! \brief QUESTION: definition?
double gcs;
//! \brief QUESTION: definition?
@@ -210,35 +212,35 @@ protected:
void allocate_vectors(C4 *c4);
public:
//! \brief QUESTION: definition?
- std::complex<double> *vh;
+ dcomplex *vh;
//! \brief QUESTION: definition?
- std::complex<double> *vj0;
+ dcomplex *vj0;
//! \brief QUESTION: definition?
- std::complex<double> *vj;
+ dcomplex *vj;
//! \brief QUESTION: definition?
- std::complex<double> *vyhj;
+ dcomplex *vyhj;
//! \brief QUESTION: definition?
- std::complex<double> *vyj0;
+ dcomplex *vyj0;
//! \brief QUESTION: definition?
- std::complex<double> **am0m;
+ dcomplex **am0m;
//! \brief QUESTION: definition?
- std::complex<double> *vint;
+ dcomplex *vint;
//! \brief QUESTION: definition?
- std::complex<double> *vintm;
+ dcomplex *vintm;
//! \brief QUESTION: definition?
- std::complex<double> **vints;
+ dcomplex **vints;
//! \brief QUESTION: definition?
- std::complex<double> *vintt;
+ dcomplex *vintt;
//! \brief QUESTION: definition?
- std::complex<double> **fsac;
+ dcomplex **fsac;
//! \brief QUESTION: definition?
- std::complex<double> **sac;
+ dcomplex **sac;
//! \brief QUESTION: definition?
- std::complex<double> **fsacm;
+ dcomplex **fsacm;
//! \brief QUESTION: definition?
double *scsc;
//! \brief QUESTION: definition?
- std::complex<double> *scscp;
+ dcomplex *scscp;
//! \brief QUESTION: definition?
double *ecsc;
//! \brief QUESTION: definition?
@@ -246,11 +248,11 @@ public:
//! \brief QUESTION: definition?
double *scscm;
//! \brief QUESTION: definition?
- std::complex<double> *ecscp;
+ dcomplex *ecscp;
//! \brief QUESTION: definition?
- std::complex<double> *scscpm;
+ dcomplex *scscpm;
//! \brief QUESTION: definition?
- std::complex<double> *ecscpm;
+ dcomplex *ecscpm;
//! \brief QUESTION: definition?
double *v3j0;
//! \brief QUESTION: definition?
@@ -296,11 +298,11 @@ protected:
int sam_size_0;
public:
//! \brief QUESTION: definition?
- std::complex<double> **gis;
+ dcomplex **gis;
//! \brief QUESTION: definition?
- std::complex<double> **gls;
+ dcomplex **gls;
//! \brief QUESTION: definition?
- std::complex<double> **sam;
+ dcomplex **sam;
/*! \brief C9 instance constructor.
*
diff --git a/src/include/Configuration.h b/src/include/Configuration.h
index f3c4df91bcc64e6673865b739b5f7dd9a06af5a3..e8e3b8951b4b01538cbf4f541b13b100d84e3e94 100644
--- a/src/include/Configuration.h
+++ b/src/include/Configuration.h
@@ -1,3 +1,5 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file Configuration.h
*
* \brief Configuration data structures.
@@ -168,7 +170,7 @@ class ScattererConfiguration {
friend void sphere(std::string, std::string, std::string);
protected:
//! \brief Matrix of dielectric parameters with size [NON_TRANS_LAYERS x N_SPHERES x N_SCALES].
- std::complex<double> ***dc0_matrix;
+ dcomplex ***dc0_matrix;
//! \brief Vector of sphere radii expressed in m, with size [N_SPHERES].
double *radii_of_spheres;
//! \brief Matrix of fractional transition radii with size [N_SPHERES x LAYERS].
@@ -258,7 +260,7 @@ public:
* \param dielectric_func_type: `int` Type of dielectric function definition (=0 for constant,
* \>0 as function of scale parameter, <0 for functions at XIP value and XI is scale factor
* for dimensions).
- * \param dc_matrix: Matrix of reference dielectric constants.
+ * \param dc_matrix: `complex double ***` Matrix of reference dielectric constants.
* \param has_external: `bool` Flag to set whether to add an external spherical layer.
* \param exdc: `double` External medium dielectric constant.
* \param wp: `double` wp
@@ -274,7 +276,7 @@ public:
int *nshl_vector,
double **rcf_vector,
int dielectric_func_type,
- std::complex<double> ***dc_matrix,
+ dcomplex ***dc_matrix,
bool has_external,
double exdc,
double wp,
diff --git a/src/include/TransitionMatrix.h b/src/include/TransitionMatrix.h
index 447457f43eb009556fcc422b75e1fd171f636ed0..3966abdab0b48463eaefb256dd6d128fe84e4e75 100644
--- a/src/include/TransitionMatrix.h
+++ b/src/include/TransitionMatrix.h
@@ -1,3 +1,5 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file TransitionMatrix.h
*
* \brief Representation of the Transition Matrix.
@@ -19,7 +21,7 @@ class TransitionMatrix {
//! External medium refractive index.
double exri;
//! Vectorized matrix elements.
- std::complex<double> *elements;
+ dcomplex *elements;
//! Sphere radius.
double sphere_radius;
//! Matrix shape
@@ -66,11 +68,11 @@ class TransitionMatrix {
* \param _lm: `int` Maximum field expansion order.
* \param _vk: `double` Wave number in scale units.
* \param _exri: `double` External medium refractive index.
- * \param _elements: `complex<double> *` Vectorized elements of the matrix.
+ * \param _elements: `complex double *` Vectorized elements of the matrix.
* \param _radius: `double` Radius for the single sphere case (defaults to 0.0).
*/
TransitionMatrix(
- int _is, int _lm, double _vk, double _exri, std::complex<double> *_elements,
+ int _is, int _lm, double _vk, double _exri, dcomplex *_elements,
double _radius=0.0
);
@@ -82,13 +84,13 @@ class TransitionMatrix {
* \param _lm: `int` Maximum field expansion order.
* \param _vk: `double` Wave number in scale units.
* \param _exri: `double` External medium refractive index.
- * \param _rmi: Matrix of complex.
- * \param _rei: Matrix of complex.
+ * \param _rmi: `complex double **`
+ * \param _rei: `complex double **`
* \param _sphere_radius: `double` Radius of the sphere.
*/
TransitionMatrix(
- int _lm, double _vk, double _exri, std::complex<double> **_rmi,
- std::complex<double> **_rei, double _sphere_radius
+ int _lm, double _vk, double _exri, dcomplex **_rmi,
+ dcomplex **_rei, double _sphere_radius
);
/*! \brief Transition Matrix instance constructor for a cluster of spheres.
@@ -100,9 +102,9 @@ class TransitionMatrix {
* \param _lm: `int` Maximum field expansion order.
* \param _vk: `double` Wave number in scale units.
* \param _exri: `double` External medium refractive index.
- * \param _am0m: Matrix of complex.
+ * \param _am0m: `complex double **`
*/
- TransitionMatrix(int _nlemt, int _lm, double _vk, double _exri, std::complex<double> **_am0m);
+ TransitionMatrix(int _nlemt, int _lm, double _vk, double _exri, dcomplex **_am0m);
/*! \brief Transition Matrix instance destroyer.
*/
diff --git a/src/include/algebraic.h b/src/include/algebraic.h
new file mode 100644
index 0000000000000000000000000000000000000000..01437ac9d3f21d49610ae6921421e5933ebefc42
--- /dev/null
+++ b/src/include/algebraic.h
@@ -0,0 +1,42 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
+/*! \file algebraic.h
+ *
+ * \brief Declaration of algebraic functions with different call-backs.
+ *
+ * In principle, the system that runs NP_TMcode may offer various types of
+ * optimized features, such as multi-core or multi-node scaling, GPU offload,
+ * or external libraries. This header collects a set of functions that can
+ * perform standard algebraic operations choosing the most optimized available
+ * system as a call-back. If no optimization is detected, eventually the
+ * legacy serial function implementation is used as a fall-back.
+ */
+
+#ifdef USE_LAPACK
+#ifdef USE_MKL
+#include <mkl_lapacke.h>
+#else
+#include <lapacke.h>
+#endif
+#else
+#define lapack_int int64_t
+#endif
+
+#ifndef INCLUDE_ALGEBRAIC_H_
+#define INCLUDE_ALGEBRAIC_H_
+
+#ifndef np_int
+#define np_int lapack_int
+#endif
+
+/*! \brief Perform in-place matrix inversion.
+ *
+ * \param mat: `complex double **` The matrix to be inverted (must be a square matrix).
+ * \param size: `np_int` The size of the matrix (i.e. the number of its rows or columns).
+ * \param ier: `int &` Reference to an integer variable for returning a result flag.
+ * \param max_size: `np_int` The maximum expected size (required by some call-backs,
+ * optional, defaults to 0).
+ */
+void invert_matrix(dcomplex **mat, np_int size, int &ier, np_int max_size=0);
+
+#endif
diff --git a/src/include/clu_subs.h b/src/include/clu_subs.h
index 49b89e20875919f916cea21a269fa7b210da3a9f..6852a46906998e878b0021dc5d606ac8500e8322 100644
--- a/src/include/clu_subs.h
+++ b/src/include/clu_subs.h
@@ -1,3 +1,5 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file clu_subs.h
*
* \brief C++ porting of CLU functions and subroutines.
@@ -15,6 +17,14 @@
#ifndef INCLUDE_CLU_SUBS_H_
#define INCLUDE_CLU_SUBS_H_
+#ifndef np_int
+#ifndef lapack_int
+#define np_int int64_t
+#else
+#define np_int lapack_int
+#endif
+#endif
+
/*! \brief Compute the asymmetry-corrected scattering cross-section of a cluster.
*
* This function computes the product between the geometrical asymmetry parameter and
@@ -23,15 +33,15 @@
*
* \param zpv: `double ****`
* \param le: `int`
- * \param am0m: Matrix of complex.
- * \param w: Matrix of complex.
+ * \param am0m: `complex double **`
+ * \param w: `complex double **`
* \param sqk: `double`
* \param gapr: `double **`
- * \param gapp: Matrix of complex.
+ * \param gapp: `complex double **`
*/
void apc(
- double ****zpv, int le, std::complex<double> **am0m, std::complex<double> **w,
- double sqk, double **gapr, std::complex<double> **gapp
+ double ****zpv, int le, dcomplex **am0m, dcomplex **w,
+ double sqk, double **gapr, dcomplex **gapp
);
/*! \brief Compute the asymmetry-corrected scattering cross-section under random average
@@ -42,32 +52,29 @@ void apc(
*
* \param zpv: `double ****`
* \param le: `int`
- * \param am0m: Matrix of complex.
+ * \param am0m: `complex double **`
* \param inpol: `int` Polarization type.
* \param sqk: `double`
* \param gaprm: `double **`
- * \param gappm: Matrix of complex.
+ * \param gappm: `complex double **`
*/
void apcra(
- double ****zpv, const int le, std::complex<double> **am0m, int inpol, double sqk,
- double **gaprm, std::complex<double> **gappm
+ double ****zpv, const int le, dcomplex **am0m, int inpol, double sqk,
+ double **gaprm, dcomplex **gappm
);
/*! \brief Complex inner product.
*
* This function performs the complex inner product. It is used by `lucin()`.
*
- * \param z: `complex<double>`
- * \param am: Matrix of complex.
+ * \param z: `complex double`
+ * \param am: `complex double **`
* \param i: `int`
* \param jf: `int`
* \param k: `int`
* \param nj: `int`
*/
-std::complex<double> cdtp(
- std::complex<double> z, std::complex<double> **am, int i, int jf,
- int k, int nj
-);
+dcomplex cdtp(dcomplex z, dcomplex **am, int i, int jf, int k, int nj);
/*! \brief C++ porting of CGEV. QUESTION: description?
*
@@ -84,13 +91,13 @@ double cgev(int ipamo, int mu, int l, int m);
* This function constructs the multi-centered M-matrix of the cluster, according
* to Eq. (5.28) of Borghese, Denti & Saija (2007).
*
- * \param am: Matrix of complex.
+ * \param am: `complex double **`
* \param c1: `C1 *`
* \param c1ao: `C1_AddOns *`
* \param c4: `C4 *`
* \param c6: `C6 *`
*/
-void cms(std::complex<double> **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6);
+void cms(dcomplex **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6);
/*! \brief Compute orientation-averaged scattered field intensity.
*
@@ -124,9 +131,9 @@ void crsm1(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6);
* \param c4: `C4 *`
* \param c6: `C6 *`
*/
-std::complex<double> ghit(
- int ihi, int ipamo, int nbl, int l1, int m1, int l2, int m2, C1 *c1,
- C1_AddOns *c1ao, C4 *c4, C6 *c6
+dcomplex ghit(
+ int ihi, int ipamo, int nbl, int l1, int m1, int l2, int m2, C1 *c1,
+ C1_AddOns *c1ao, C4 *c4, C6 *c6
);
/*! \brief Compute Hankel funtion and Bessel functions.
@@ -144,7 +151,7 @@ std::complex<double> ghit(
* \param c4: `C4 *`
*/
void hjv(
- double exri, double vk, int &jer, int &lcalc, std::complex<double> &arg,
+ double exri, double vk, int &jer, int &lcalc, dcomplex &arg,
C1 *c1, C1_AddOns *c1ao, C4 *c4
);
@@ -153,12 +160,12 @@ void hjv(
* This function performs the inversion of the multi-centered M-matrix through
* LU decomposition. See Eq. (5.29) in Borghese, Denti & Saija (2007).
*
- * \param am: Matrix of complex.
- * \param nddmst: `const int`
- * \param n: `int`
+ * \param am: `complex double **`
+ * \param nddmst: `const int64_t`
+ * \param n: `int64_t`
* \param ier: `int &`
*/
-void lucin(std::complex<double> **am, const int nddmst, int n, int &ier);
+void lucin(dcomplex **am, const np_int nddmst, np_int n, int &ier);
/*! \brief Compute the average extinction cross-section.
*
@@ -172,7 +179,7 @@ void lucin(std::complex<double> **am, const int nddmst, int n, int &ier);
* \param cextlr: `double **`
* \param cext: `double **`
*/
-void mextc(double vk, double exri, std::complex<double> **fsac, double **cextlr, double **cext);
+void mextc(double vk, double exri, dcomplex **fsac, double **cextlr, double **cext);
/*! \brief Compute cross-sections and forward scattering amplitude for the cluster.
*
@@ -266,16 +273,16 @@ void r3jmr(int j1, int j2, int j3, int m1, C6 *c6);
* in Borghese, Denti & Saija (2007).
*
* \param le: `int`
- * \param am0m: Matrix of complex.
- * \param w: Matrix of complex.
+ * \param am0m: `complex double **`
+ * \param w: `complex double **`
* \param tqce: `double **`
- * \param tqcpe: Matrix of complex.
+ * \param tqcpe: `complex double **`
* \param tqcs: `double **`
- * \param tqcps: Matrix of complex.
+ * \param tqcps: `complex double **`
*/
void raba(
- int le, std::complex<double> **am0m, std::complex<double> **w, double **tqce,
- std::complex<double> **tqcpe, double **tqcs, std::complex<double> **tqcps
+ int le, dcomplex **am0m, dcomplex **w, double **tqce,
+ dcomplex **tqcpe, double **tqcs, dcomplex **tqcps
);
/*! \brief Compute the radiation force Cartesian components.
@@ -383,13 +390,13 @@ void tqr(
* This function computes the single-centered inverrted M-matrix appearing in Eq. (5.28)
* of Borghese, Denti & Saija (2007).
*
- * \param am: Matrix of complex.
+ * \param am: `complex double **`
* \param c1: `C1 *` Pointer to a C1 instance.
* \param c1ao: `C1_AddOns *` Pointer to C1_AddOns instance.
* \param c4: `C4 *` Pointer to a C4 structure.
* \param c6: `C6 *` Pointer to a C6 instance.
* \param c9: `C9 *` Pointer to a C9 instance.
*/
-void ztm(std::complex<double> **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6, C9 * c9);
+void ztm(dcomplex **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6, C9 * c9);
#endif
diff --git a/src/include/lapack_calls.h b/src/include/lapack_calls.h
new file mode 100644
index 0000000000000000000000000000000000000000..6f2f365bbc027bf2c4037ccbc2dcbe33403874bf
--- /dev/null
+++ b/src/include/lapack_calls.h
@@ -0,0 +1,23 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
+/*! \file lapack_calss.h
+ *
+ * \brief C++ interface to LAPACK calls.
+ *
+ */
+
+#ifndef INCLUDE_LAPACK_CALLS_H_
+#define INCLUDE_LAPACK_CALLS_H_
+
+/*! \brief Invert a complex matrix with double precision elements.
+ *
+ * Use LAPACKE64 to perform an in-place matrix inversion for a complex
+ * matrix with double precision elements.
+ *
+ * \param mat: Matrix of complex. The matrix to be inverted.
+ * \param n: `lapack_int` The number of rows and columns of the [n x n] matrix.
+ * \param jer: `int &` Reference to an integer return flag.
+ */
+void zinvert(dcomplex **mat, lapack_int n, int &jer);
+
+#endif
diff --git a/src/include/sph_subs.h b/src/include/sph_subs.h
index 2aafef2233caa43a8f76377afc65db380f056c45..427ae3a63af00fa92a1d2df86a8babcbaa91746c 100644
--- a/src/include/sph_subs.h
+++ b/src/include/sph_subs.h
@@ -1,3 +1,5 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file sph_subs.h
*
* \brief C++ porting of SPH functions and subroutines.
@@ -35,11 +37,11 @@ void aps(double ****zpv, int li, int nsph, C1 *c1, double sqk, double *gaps);
* backward recurrence. This is the `CSPHJ` implementation of the `specfun` library.
*
* \param n: `int` Order of the function.
- * \param z: `complex<double>` Argument of the function.
+ * \param z: `complex double` Argument of the function.
* \param nm: `int &` Highest computed order.
- * \param csj: Vector of complex. The desired function \f$j\f$.
+ * \param csj: `complex double *` The desired function \f$j\f$.
*/
-void cbf(int n, std::complex<double> z, int &nm, std::complex<double> csj[]);
+void cbf(int n, dcomplex z, int &nm, dcomplex *csj);
/*! \brief Clebsch-Gordan coefficients.
*
@@ -56,10 +58,10 @@ double cg1(int lmpml, int mu, int l, int m);
/*! \brief Conjugate of a double precision complex number
*
- * \param z: `complex<double>` The input complex number.
- * \return result: `complex<double>` The conjugate of the input number.
+ * \param z: `complex double` The input complex number.
+ * \return result: `complex double` The conjugate of the input number.
*/
-std::complex<double> dconjg(std::complex<double> z);
+dcomplex dconjg(dcomplex z);
/*! \brief Comute the continuous variation of the refractive index and of its radial derivative.
*
@@ -95,11 +97,11 @@ void diel(int npntmo, int ns, int i, int ic, double vk, C1 *c1, C2 *c2);
* \param c2: `C2 *` Pointer to a `C2` data structure.
* \param jer: `int &` Reference to integer error code variable.
* \param lcalc: `int &` Reference to integer variable recording the maximum expansion order accounted for.
- * \param arg: `complex<double> &`
+ * \param arg: `complex double &`
*/
void dme(
int li, int i, int npnt, int npntts, double vk, double exdc, double exri,
- C1 *c1, C2 *c2, int &jer, int &lcalc, std::complex<double> &arg
+ C1 *c1, C2 *c2, int &jer, int &lcalc, dcomplex &arg
);
/*! \brief Bessel function calculation control parameters.
@@ -118,11 +120,11 @@ double envj(int n, double x);
* This function computes the Mueller Transformation Matrix, or Phase Matrix. See
* Sec. 2.8.1 of Borghese, Denti & Saija (2007).
*
- * \param vint: Vector of complex.
+ * \param vint: `complex double *`
* \param cmullr: `double **`
* \param cmul: `double **`
*/
-void mmulc(std::complex<double> *vint, double **cmullr, double **cmul);
+void mmulc(dcomplex *vint, double **cmullr, double **cmul);
/*! \brief Starting point for Bessel function magnitude.
*
@@ -168,14 +170,14 @@ void orunve(double *u1, double *u2, double *u3, int iorth, double torth);
*
* \param up: `double *`
* \param un: `double *`
- * \param ylm: Vector of complex. Field polar spherical harmonics.
+ * \param ylm: `complex double *` Field polar spherical harmonics.
* \param inpol: `int` Incident field polarization type (0 - linear; 1 - circular).
* \param lw: `int`
* \param isq: `int`
* \param c1: `C1 *` Pointer to a `C1` data structure.
*/
void pwma(
- double *up, double *un, std::complex<double> *ylm, int inpol, int lw,
+ double *up, double *un, dcomplex *ylm, int inpol, int lw,
int isq, C1 *c1
);
@@ -189,14 +191,14 @@ void pwma(
* \param li: `int` Maximum field expansion order.
* \param nsph: `int` Number of spheres.
* \param c1: `C1 *` Pointer to `C1` data structure.
- * \param tqse: Matrix of complex.
- * \param tqspe: Matrix of complex.
- * \param tqss: Matrix of complex.
- * \param tqsps: Matrix of complex.
+ * \param tqse: `double **`
+ * \param tqspe: `complex double **`
+ * \param tqss: `double **`
+ * \param tqsps: `complex double **`
*/
void rabas(
- int inpol, int li, int nsph, C1 *c1, double **tqse, std::complex<double> **tqspe,
- double **tqss, std::complex<double> **tqsps
+ int inpol, int li, int nsph, C1 *c1, double **tqse, dcomplex **tqspe,
+ double **tqss, dcomplex **tqsps
);
/*! \brief Real Bessel Function.
@@ -219,18 +221,17 @@ void rbf(int n, double x, int &nm, double sj[]);
*
* \param npntmo: `int`
* \param step: `double`
- * \param dcc: `complex<double>`
+ * \param dcc: `complex double`
* \param x: `double &`
* \param lpo: `int`
- * \param y1: `complex<double> &`
- * \param y2: `complex<double> &`
- * \param dy1: `complex<double> &`
- * \param dy2: `complex<double> &`
+ * \param y1: `complex double &`
+ * \param y2: `complex double &`
+ * \param dy1: `complex double &`
+ * \param dy2: `complex double &`
*/
void rkc(
- int npntmo, double step, std::complex<double> dcc, double &x, int lpo,
- std::complex<double> &y1, std::complex<double> &y2, std::complex<double> &dy1,
- std::complex<double> &dy2
+ int npntmo, double step, dcomplex dcc, double &x, int lpo,
+ dcomplex &y1, dcomplex &y2, dcomplex &dy1, dcomplex &dy2
);
/*! \brief Transition layer radial function and derivative.
@@ -242,16 +243,15 @@ void rkc(
* \param step: `double`
* \param x: `double &`
* \param lpo: `int`
- * \param y1: `complex<double> &`
- * \param y2: `complex<double> &`
- * \param dy1: `complex<double> &`
- * \param dy2: `complex<double> &`
+ * \param y1: `complex double &`
+ * \param y2: `complex double &`
+ * \param dy1: `complex double &`
+ * \param dy2: `complex double &`
* \param c2: `C2 *` Pointer to a `C2` data structure.
*/
void rkt(
- int npntmo, double step, double &x, int lpo, std::complex<double> &y1,
- std::complex<double> &y2, std::complex<double> &dy1, std::complex<double> &dy2,
- C2 *c2
+ int npntmo, double step, double &x, int lpo, dcomplex &y1,
+ dcomplex &y2, dcomplex &dy1, dcomplex &dy2, C2 *c2
);
/*! \brief Spherical Bessel functions.
@@ -262,9 +262,9 @@ void rkt(
* \param n: `int` Order of the function (from 0 up).
* \param x: `double` Argumento of the function (\f$x > 0\f$).
* \param nm: `int &` Highest computed order.
- * \param sy: `double[]` The desired function \f$y\f$.
+ * \param sy: `double *` The desired function \f$y\f$.
*/
-void rnf(int n, double x, int &nm, double sy[]);
+void rnf(int n, double x, int &nm, double *sy);
/*! \brief Spherical harmonics for given direction.
*
@@ -276,11 +276,11 @@ void rnf(int n, double x, int &nm, double sy[]);
* \param cosrph: `double` Cosine of direction's azimuth.
* \param sinrph: `double` Sine of direction's azimuth.
* \param ll: `int` L value expansion order.
- * \param ylm: Vector of complex. The requested spherical harmonics.
+ * \param ylm: `complex double *` The requested spherical harmonics.
*/
void sphar(
double cosrth, double sinrth, double cosrph, double sinrph,
- int ll, std::complex<double> *ylm
+ int ll, dcomplex *ylm
);
/*! \brief Compute scattering, absorption and extinction cross-sections.
@@ -288,14 +288,14 @@ void sphar(
* This function computes the scattering, absorption and extinction cross-sections in terms
* of Forward Scattering Amplitudes. See Sec. 4.2.1 in Borghese, Denti & Saija (2007).
*
- * \param tfsas: `complex<double> &`
+ * \param tfsas: `complex double &`
* \param nsph: `int` Number of spheres.
* \param lm: `int` Maximum field expansion order.
* \param vk: `double` Wave number in scale units.
* \param exri: `double` External medium refractive index.
* \param c1: `C1 *` Pointer to a `C1` data structure.
*/
-void sscr0(std::complex<double> &tfsas, int nsph, int lm, double vk, double exri, C1 *c1);
+void sscr0(dcomplex &tfsas, int nsph, int lm, double vk, double exri, C1 *c1);
/*! \brief C++ Compute the scattering amplitude and the scattered field intensity.
*
diff --git a/src/include/tfrfme.h b/src/include/tfrfme.h
index ba2954a742b5458f8cabb1ca64729befef232e20..2e55c4822b1d186bb7a8094e6ba89aafa46f176c 100644
--- a/src/include/tfrfme.h
+++ b/src/include/tfrfme.h
@@ -1,3 +1,5 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file tfrfme.h
*
* \brief Representation of the trapping calculation objects.
@@ -18,7 +20,7 @@ protected:
int nlmmt;
//! QUESTION: definition?
- std::complex<double> *wk;
+ dcomplex *wk;
/*! \brief Load a Swap1 instance from a HDF5 binary file.
*
@@ -60,9 +62,9 @@ public:
/*! \brief Append an element at the end of the vector.
*
- * \param value: `complex<double>` The value to be added to the vector.
+ * \param value: `complex double` The value to be added to the vector.
*/
- void append(std::complex<double> value) { wk[last_index++] = value; }
+ void append(dcomplex value) { wk[last_index++] = value; }
/*! \brief Load a Swap1 instance from binary file.
*
@@ -83,9 +85,9 @@ public:
/*! \brief Get the pointer to the WK vector.
*
- * \return value: `complex<double> *` Memory address of the WK vector.
+ * \return value: `complex double *` Memory address of the WK vector.
*/
- std::complex<double> *get_vector() { return wk; }
+ dcomplex *get_vector() { return wk; }
/*! \brief Bring the pointer to the next element at the start of vector.
*/
@@ -311,7 +313,7 @@ protected:
//! Vector of computed z positions
double *zv;
//! QUESTION: definition?
- std::complex<double> **wsum;
+ dcomplex **wsum;
/*! \brief Load a configuration instance from a HDF5 binary file.
*
@@ -371,9 +373,9 @@ public:
/*! \brief Get the pointer to the WSUM matrix.
*
- * \return value: `complex<double> **` Pointer to the WSUM matrix.
+ * \return value: `complex double **` Pointer to the WSUM matrix.
*/
- std::complex<double> **get_matrix() { return wsum; }
+ dcomplex **get_matrix() { return wsum; }
/*! \brief Calculate the necessary amount of memory to create a new instance.
*
diff --git a/src/include/tra_subs.h b/src/include/tra_subs.h
index c32aa9ae9bffb58387d936afcf6547c12713be21..37e728507fe71f9444d7f1fd405044b2aad6e3e4 100644
--- a/src/include/tra_subs.h
+++ b/src/include/tra_subs.h
@@ -1,3 +1,5 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file tra_subs.h
*
* \brief C++ porting of TRAPPING functions and subroutines.
@@ -48,17 +50,14 @@ struct CCR {
/*! C++ porting of CAMP
*
- * \param ac: Vector of complex. QUESTION: definition?
- * \param am0m: Matrix of complex. QUESTION: definition?
- * \param ws: Vector of complex. QUESTION: definition?
+ * \param ac: `complex double *` QUESTION: definition?
+ * \param am0m: `complex double **` QUESTION: definition?
+ * \param ws: `complex double *` QUESTION: definition?
* \param cil: `CIL *` Pointer to a CIL structure.
*
* This function builds the AC vector using AM0M and WS.
*/
-void camp(
- std::complex<double> *ac, std::complex<double> **am0m, std::complex<double> *ws,
- CIL *cil
-);
+void camp(dcomplex *ac, dcomplex **am0m, dcomplex *ws, CIL *cil);
/*! C++ porting of CZAMP
*
@@ -70,26 +69,33 @@ void camp(
*
* This function builds the AC vector using AMD, INDAM and WS.
*/
-void czamp(
- std::complex<double> *ac, std::complex<double> **amd, int **indam,
- std::complex<double> *ws, CIL *cil
-);
+void czamp(dcomplex *ac, dcomplex **amd, int **indam, dcomplex *ws, CIL *cil);
/*! C++ porting of FFRF
*
* \param zpv: `double ****`. QUESTION: definition?
- * \param ac: Vector of complex. QUESTION: definition?
- * \param ws: Vector of complex. QUESTION: definition?
+ * \param ac: `complex double *` QUESTION: definition?
+ * \param ws: `complex double *` QUESTION: definition?
* \param fffe: `double *`. QUESTION: definition?
* \param fffs: `double *`. QUESTION: definition?
* \param cil: `CIL *` Pointer to a CIL structure.
* \param ccr: `CCR *` Pointer to a CCR structure.
*/
void ffrf(
- double ****zpv, std::complex<double> *ac, std::complex<double> *ws, double *fffe,
+ double ****zpv, dcomplex *ac, dcomplex *ws, double *fffe,
double *fffs, CIL *cil, CCR *ccr
);
+/*! C++ porting of FFRT
+ *
+ * \param ac: `complex double *` QUESTION: definition?
+ * \param ws: `complex double *` QUESTION: definition?
+ * \param ffte: `double *`. QUESTION: definition?
+ * \param ffts: `double *`. QUESTION: definition?
+ * \param cil: `CIL *` Pointer to a CIL structure.
+ */
+void ffrt(dcomplex *ac, dcomplex *ws, double *ffte, double *ffts, CIL *cil);
+
/*! C++ porting of FRFMER
*
* \param nkv: `int` QUESTION: definition?
@@ -105,68 +111,50 @@ void ffrf(
* \param pmf: `double` QUESTION: definition?
* \param tt1: `Swap1 *` Pointer to first swap object.
* \param tt2: `Swap2 *` Pointer to second swap object.
+ * \return wk: `complex double *` QUESTION: definition?
*/
-std::complex<double> *frfmer(
+dcomplex *frfmer(
int nkv, double vkm, double vknmx, double apfafa, double tra,
double spd, double rir, double ftcn, int le, int lmode, double pmf,
Swap1 *tt1, Swap2 *tt2
);
-/*! C++ porting of FFRT
- *
- * \param ac: Vector of complex. QUESTION: definition?
- * \param ws: Vector of complex. QUESTION: definition?
- * \param ffte: `double *`. QUESTION: definition?
- * \param ffts: `double *`. QUESTION: definition?
- * \param cil: `CIL *` Pointer to a CIL structure.
- */
-void ffrt(
- std::complex<double> *ac, std::complex<double> *ws, double *ffte, double *ffts,
- CIL *cil
-);
-
/*! C++ porting of PWMALP
*
- * \param w: Matrix of complex. QUESTION: definition?
+ * \param w: `complex double *` QUESTION: definition?
* \param up: `double *`
* \param un: `double *`
- * \param ylm: Vector of complex
+ * \param ylm: `complex double *` Field vector spherical harmonics.
* \param lw: `int`
*/
-void pwmalp(std::complex<double> **w, double *up, double *un, std::complex<double> *ylm, int lw);
+void pwmalp(dcomplex **w, double *up, double *un, dcomplex *ylm, int lw);
/*! C++ porting of SAMP
*
- * \param ac: Vector of complex. QUESTION: definition?
- * \param tmsm: Vector of complex. QUESTION: definition?
- * \param tmse: Vector of complex. QUESTION: definition?
- * \param ws: Vector of complex. QUESTION: definition?
+ * \param ac: `complex double *` QUESTION: definition?
+ * \param tmsm: `complex double *` QUESTION: definition?
+ * \param tmse: `complex double *` QUESTION: definition?
+ * \param ws: `complex double *` QUESTION: definition?
* \param cil: `CIL *` Pointer to a CIL structure.
*
* This function builds the AC vector using TMSM, TMSE and WS.
*/
-void samp(
- std::complex<double> *ac, std::complex<double> *tmsm, std::complex<double> *tmse,
- std::complex<double> *ws, CIL *cil
-);
+void samp(dcomplex *ac, dcomplex *tmsm, dcomplex *tmse, dcomplex *ws, CIL *cil);
/*! C++ porting of SAMPOA
*
- * \param ac: Vector of complex. QUESTION: definition?
- * \param tms: Matrix of complex. QUESTION: definition?
- * \param ws: Vector of complex. QUESTION: definition?
+ * \param ac: `complex double *` QUESTION: definition?
+ * \param tms: `complex double **` QUESTION: definition?
+ * \param ws: `complex double *` QUESTION: definition?
* \param cil: `CIL *` Pointer to a CIL structure.
*
* This function builds the AC vector using TMS and WS.
*/
-void sampoa(
- std::complex<double> *ac, std::complex<double> **tms, std::complex<double> *ws,
- CIL *cil
-);
+void sampoa(dcomplex *ac, dcomplex **tms, dcomplex *ws, CIL *cil);
/*! C++ porting of WAMFF
*
- * \param wk: `complex<double> *` QUESTION: definition?
+ * \param wk: `complex double *` QUESTION: definition?
* \param x: `double`
* \param y: `double`
* \param z: `double`
@@ -180,6 +168,6 @@ void sampoa(
* \param pmf: `double` QUESTION: definition?
*/
void wamff(
- std::complex<double> *wk, double x, double y, double z, int lm, double apfafa,
+ dcomplex *wk, double x, double y, double z, int lm, double apfafa,
double tra, double spd, double rir, double ftcn, int lmode, double pmf
);
diff --git a/src/include/types.h b/src/include/types.h
new file mode 100644
index 0000000000000000000000000000000000000000..0271769a73600adc8de3ff6dfe1cf4d9d0f8a49f
--- /dev/null
+++ b/src/include/types.h
@@ -0,0 +1,28 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
+/*! \file types.h
+ *
+ * \brief Definition of fundamental types in use.
+ */
+
+#ifndef INCLUDE_TYPES_H_
+#define INCLUDE_TYPES_H_
+
+#include <complex.h>
+
+typedef __complex__ double dcomplex;
+
+/*! \brief Get the real part of a complex number.
+ *
+ * \param z: `complex double` The argument of the function.
+ * \return rz: `double` The real part of the argument.
+ */
+double real(dcomplex z);
+
+/*! \brief Get the imaginary part of a complex number.
+ *
+ * \param z: `complex double` The argument of the function.
+ * \return iz: `double` The imaginary part of the argument.
+ */
+double imag(dcomplex z);
+#endif
diff --git a/src/libnptm/Commons.cpp b/src/libnptm/Commons.cpp
index d19e62b58843443c3c1540d43f353b05e55bd025..fc9498a839ede48183ce665b3fd15c01d8e404cb 100644
--- a/src/libnptm/Commons.cpp
+++ b/src/libnptm/Commons.cpp
@@ -1,3 +1,5 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file Commons.cpp
*
* DEVELOPMENT NOTE:
@@ -10,38 +12,38 @@
* to the configuration objects. These, on their turn, need to
* expose methods to access the relevant data in read-only mode.
*/
-#include <complex>
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
#ifndef INCLUDE_COMMONS_H
#include "../include/Commons.h"
#endif
-using namespace std;
-
C1::C1(int ns, int l_max, int *_nshl, int *_iog) {
nlmmt = 2 * (l_max * (l_max + 2));
nsph = ns;
lm = l_max;
- rmi = new complex<double>*[lm];
- rei = new complex<double>*[lm];
+ rmi = new dcomplex*[lm];
+ rei = new dcomplex*[lm];
for (int ri = 0; ri < lm; ri++) {
- rmi[ri] = new complex<double>[nsph]();
- rei[ri] = new complex<double>[nsph]();
+ rmi[ri] = new dcomplex[nsph]();
+ rei[ri] = new dcomplex[nsph]();
}
- w = new complex<double>*[nlmmt];
- for (int wi = 0; wi < nlmmt; wi++) w[wi] = new complex<double>[4]();
+ w = new dcomplex*[nlmmt];
+ for (int wi = 0; wi < nlmmt; wi++) w[wi] = new dcomplex[4]();
int configurations = 0;
for (int ci = 1; ci <= nsph; ci++) {
if (_iog[ci - 1] >= ci) configurations++;
}
- vints = new complex<double>*[nsph];
+ vints = new dcomplex*[nsph];
rc = new double*[configurations];
nshl = new int[configurations]();
iog = new int[nsph]();
int conf_index = 0;
for (int vi = 0; vi < nsph; vi++) {
- vints[vi] = new complex<double>[16]();
+ vints[vi] = new dcomplex[16]();
iog[vi] = _iog[vi];
if (iog[vi] >= vi + 1) {
nshl[conf_index] = _nshl[conf_index];
@@ -49,7 +51,7 @@ C1::C1(int ns, int l_max, int *_nshl, int *_iog) {
conf_index++;
}
}
- fsas = new complex<double>[nsph]();
+ fsas = new dcomplex[nsph]();
sscs = new double[nsph]();
sexs = new double[nsph]();
sabs = new double[nsph]();
@@ -62,11 +64,11 @@ C1::C1(int ns, int l_max, int *_nshl, int *_iog) {
rzz = new double[nsph]();
ros = new double[nsph]();
- sas = new complex<double>**[nsph];
+ sas = new dcomplex**[nsph];
for (int si = 0; si < nsph; si++) {
- sas[si] = new complex<double>*[2];
- sas[si][0] = new complex<double>[2]();
- sas[si][1] = new complex<double>[2]();
+ sas[si] = new dcomplex*[2];
+ sas[si][0] = new dcomplex[2]();
+ sas[si][1] = new dcomplex[2]();
}
}
@@ -112,32 +114,32 @@ C1_AddOns::C1_AddOns(C4 *c4) {
nsph = c4->nsph;
lmpo = c4->lmpo;
nlemt = 2 * c4->nlem;
- vh = new complex<double>[(nsph * nsph - 1) * c4->litpo]();
- vj0 = new complex<double>[nsph * c4->lmtpo]();
- vj = new complex<double>[1](); // QUESTION: is 1 really enough for a general case?
- vyhj = new complex<double>[(nsph * nsph - 1) * c4->litpos]();
- vyj0 = new complex<double>[nsph * c4->lmtpos]();
- am0m = new complex<double>*[nlemt];
+ vh = new dcomplex[(nsph * nsph - 1) * c4->litpo]();
+ vj0 = new dcomplex[nsph * c4->lmtpo]();
+ vj = new dcomplex[1](); // QUESTION: is 1 really enough for a general case?
+ vyhj = new dcomplex[(nsph * nsph - 1) * c4->litpos]();
+ vyj0 = new dcomplex[nsph * c4->lmtpos]();
+ am0m = new dcomplex*[nlemt];
for (int ai = 0; ai < nlemt; ai++) {
- am0m[ai] = new complex<double>[nlemt]();
+ am0m[ai] = new dcomplex[nlemt]();
}
- vint = new complex<double>[16](); // QUESTION: is dimension 16 generally fixed?
- vintm = new complex<double>[16]();
- vintt = new complex<double>[16]();
- vints = new complex<double>*[nsph];
- for (int vi = 0; vi < nsph; vi++) vints[vi] = new complex<double>[16]();
- fsac = new complex<double>*[2];
- sac = new complex<double>*[2];
- fsacm = new complex<double>*[2];
+ vint = new dcomplex[16](); // QUESTION: is dimension 16 generally fixed?
+ vintm = new dcomplex[16]();
+ vintt = new dcomplex[16]();
+ vints = new dcomplex*[nsph];
+ for (int vi = 0; vi < nsph; vi++) vints[vi] = new dcomplex[16]();
+ fsac = new dcomplex*[2];
+ sac = new dcomplex*[2];
+ fsacm = new dcomplex*[2];
for (int fi = 0; fi < 2; fi++) {
- fsac[fi] = new complex<double>[2]();
- sac[fi] = new complex<double>[2]();
- fsacm[fi] = new complex<double>[2]();
+ fsac[fi] = new dcomplex[2]();
+ sac[fi] = new dcomplex[2]();
+ fsacm[fi] = new dcomplex[2]();
}
- scscp = new complex<double>[2]();
- ecscp = new complex<double>[2]();
- scscpm = new complex<double>[2]();
- ecscpm = new complex<double>[2]();
+ scscp = new dcomplex[2]();
+ ecscp = new dcomplex[2]();
+ scscpm = new dcomplex[2]();
+ ecscpm = new dcomplex[2]();
allocate_vectors(c4);
sscs = new double[nsph]();
ecscm = new double[2]();
@@ -191,20 +193,20 @@ void C1_AddOns::allocate_vectors(C4 *c4) {
for (int ii = 0; ii < lmpo; ii++) ind3j[ii] = new int[c4->lm]();
// This calculates the size of vyhj
int ivy = (nsph * nsph - 1) * c4->litpos;
- vyhj = new complex<double>[ivy]();
+ vyhj = new dcomplex[ivy]();
// This calculates the size of vyj0
ivy = c4->lmtpos * c4->nsph;
- vyj0 = new complex<double>[ivy]();
+ vyj0 = new dcomplex[ivy]();
}
C2::C2(int ns, int nl, int npnt, int npntts) {
nsph = ns;
int max_n = (npnt > npntts) ? npnt : npntts;
nhspo = 2 * max_n - 1;
- ris = new complex<double>[nhspo]();
- dlri = new complex<double>[nhspo]();
- vkt = new complex<double>[nsph]();
- dc0 = new complex<double>[nl]();
+ ris = new dcomplex[nhspo]();
+ dlri = new dcomplex[nhspo]();
+ vkt = new dcomplex[nsph]();
+ dc0 = new dcomplex[nl]();
vsz = new double[nsph]();
}
@@ -217,10 +219,10 @@ C2::~C2() {
}
C3::C3() {
- tsas = new complex<double>*[2];
- tsas[0] = new complex<double>[2];
- tsas[1] = new complex<double>[2];
- tfsas = complex<double>(0.0, 0.0);
+ tsas = new dcomplex*[2];
+ tsas[0] = new dcomplex[2];
+ tsas[1] = new dcomplex[2];
+ tfsas = 0.0 + 0.0 * I;
gcs = 0.0;
scs = 0.0;
ecs = 0.0;
@@ -242,14 +244,14 @@ C6::~C6() {
}
C9::C9(int ndi, int nlem, int ndit, int nlemt) {
- gis = new complex<double>*[ndi];
- gls = new complex<double>*[ndi];
+ gis = new dcomplex*[ndi];
+ gls = new dcomplex*[ndi];
for (int gi = 0; gi < ndi; gi++) {
- gis[gi] = new complex<double>[nlem]();
- gls[gi] = new complex<double>[nlem]();
+ gis[gi] = new dcomplex[nlem]();
+ gls[gi] = new dcomplex[nlem]();
}
- sam = new complex<double>*[ndit];
- for (int si = 0; si < ndit; si++) sam[si] = new complex<double>[nlemt]();
+ sam = new dcomplex*[ndit];
+ for (int si = 0; si < ndit; si++) sam[si] = new dcomplex[nlemt]();
gis_size_0 = ndi;
sam_size_0 = ndit;
}
diff --git a/src/libnptm/Configuration.cpp b/src/libnptm/Configuration.cpp
index 04cdb5c03a0930affa172afb46d7db960166914b..1bf78b7db3e571858bee6ad893f0c575dbb33644 100644
--- a/src/libnptm/Configuration.cpp
+++ b/src/libnptm/Configuration.cpp
@@ -1,10 +1,11 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file Configuration.cpp
*
* \brief Implementation of the configuration classes.
*/
#include <cmath>
-#include <complex>
#include <cstdio>
#include <exception>
#include <fstream>
@@ -12,6 +13,10 @@
#include <regex>
#include <string>
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
+
#ifndef INCLUDE_ERRORS_H_
#include "../include/errors.h"
#endif
@@ -202,7 +207,7 @@ ScattererConfiguration::ScattererConfiguration(
int *nshl_vector,
double **rcf_vector,
int dielectric_func_type,
- complex<double> ***dc_matrix,
+ dcomplex ***dc_matrix,
bool is_external,
double ex,
double w,
@@ -477,12 +482,12 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(string dedfb_file_nam
} // ns112 loop
} // i113 loop
- complex<double> ***dc0m = new complex<double>**[configurations];
- complex<double> *dc0 = new complex<double>[configurations]();
+ dcomplex ***dc0m = new dcomplex**[configurations];
+ dcomplex *dc0 = new dcomplex[configurations]();
int dim3 = (_idfc == 0) ? nxi : 1;
for (int di = 0; di < configurations; di++) {
- dc0m[di] = new complex<double>*[nsph];
- for (int dj = 0; dj < nsph; dj++) dc0m[di][dj] = new complex<double>[dim3]();
+ dc0m[di] = new dcomplex*[nsph];
+ for (int dj = 0; dj < nsph; dj++) dc0m[di][dj] = new dcomplex[dim3]();
} // di loop
for (int jxi468 = 1; jxi468 <= nxi; jxi468++) {
if (_idfc != 0 && jxi468 > 1) continue; // jxi468 loop
@@ -504,7 +509,7 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(string dedfb_file_nam
str_number = regex_replace(str_number, regex("D"), "e");
str_number = regex_replace(str_number, regex("d"), "e");
double ival = stod(str_number);
- dc0[ic157] = complex<double>(rval, ival);
+ dc0[ic157] = rval + ival * I;
dc0m[ic157][i162 - 1][jxi468 - 1] = dc0[ic157];
} // ic157 loop
}
@@ -547,7 +552,7 @@ ScattererConfiguration* ScattererConfiguration::from_hdf5(string file_name) {
double *xi_vec;
double *ros_vector;
double **rcf_vector;
- complex<double> ***dc0m;
+ dcomplex ***dc0m;
status = hdf_file->read("NSPH", "INT32_(1)", &nsph);
status = hdf_file->read("IES", "INT32_(1)", &ies);
status = hdf_file->read("EXDC", "FLOAT64_(1)", &_exdc);
@@ -591,17 +596,17 @@ ScattererConfiguration* ScattererConfiguration::from_hdf5(string file_name) {
str_name = "DC0M";
str_type = "FLOAT64_(" + to_string(element_size) + ")";
status = hdf_file->read(str_name, str_type, elements);
- dc0m = new complex<double>**[configuration_count];
+ dc0m = new dcomplex**[configuration_count];
int dc_index = 0;
for (int di = 0; di < configuration_count; di++) {
- dc0m[di] = new complex<double>*[nsph];
+ dc0m[di] = new dcomplex*[nsph];
for (int dj = 0; dj < nsph; dj++) {
- dc0m[di][dj] = new complex<double>[dim3]();
+ dc0m[di][dj] = new dcomplex[dim3]();
for (int dk = 0; dk < dim3; dk++) {
double rval = elements[2 * dc_index];
double ival = elements[2 * dc_index + 1];
dc_index++;
- dc0m[di][dj][dk] = complex<double>(rval, ival);
+ dc0m[di][dj][dk] = rval + ival * I;
}
} // dj loop
} // di loop
@@ -638,7 +643,7 @@ ScattererConfiguration* ScattererConfiguration::from_legacy(string file_name) {
double *_xi_vec;
double *_ros_vec;
double **_rcf_vec;
- complex<double> ***_dc0m;
+ dcomplex ***_dc0m;
fstream input;
input.open(file_name.c_str(), ios::in | ios::binary);
@@ -672,11 +677,11 @@ ScattererConfiguration* ScattererConfiguration::from_legacy(string file_name) {
for (int nsi = 0; nsi < nsh; nsi++)
input.read(reinterpret_cast<char *>(&(_rcf_vec[i115 - 1][nsi])), sizeof(double));
}
- _dc0m = new complex<double>**[configurations];
+ _dc0m = new dcomplex**[configurations];
int dim3 = (_idfc == 0) ? _nxi : 1;
for (int di = 0; di < configurations; di++) {
- _dc0m[di] = new complex<double>*[_nsph];
- for (int dj = 0; dj < _nsph; dj++) _dc0m[di][dj] = new complex<double>[dim3]();
+ _dc0m[di] = new dcomplex*[_nsph];
+ for (int dj = 0; dj < _nsph; dj++) _dc0m[di][dj] = new dcomplex[dim3]();
} // di loop
for (int jxi468 = 1; jxi468 <= _nxi; jxi468++) {
if (_idfc != 0 && jxi468 > 1) continue;
@@ -689,7 +694,7 @@ ScattererConfiguration* ScattererConfiguration::from_legacy(string file_name) {
double rval, ival;
input.read(reinterpret_cast<char *>(&rval), sizeof(double));
input.read(reinterpret_cast<char *>(&ival), sizeof(double));
- _dc0m[i157][i162 - 1][jxi468 - 1] = complex<double>(rval, ival);
+ _dc0m[i157][i162 - 1][jxi468 - 1] = rval + ival * I;
}
}
}
@@ -771,7 +776,7 @@ void ScattererConfiguration::print() {
printf(" [");
for (int k = 0; k < number_of_scales; k++) {
if (idfc != 0 and k > 0) continue;
- printf("\t%lg + i(%lg)", dc0_matrix[i][j][k].real(), dc0_matrix[i][j][k].imag());
+ printf("\t%lg + i(%lg)", real(dc0_matrix[i][j][k]), imag(dc0_matrix[i][j][k]));
}
printf("\t]\n");
}
@@ -862,8 +867,8 @@ void ScattererConfiguration::write_hdf5(string file_name) {
if (i162 == 1) ici = ici + ies;
for (int i157 = 0; i157 < ici; i157++) {
double dc0_real, dc0_imag;
- dc0_real = dc0_matrix[i157][i162 - 1][jxi468 - 1].real();
- dc0_imag = dc0_matrix[i157][i162 - 1][jxi468 - 1].imag();
+ dc0_real = real(dc0_matrix[i157][i162 - 1][jxi468 - 1]);
+ dc0_imag = imag(dc0_matrix[i157][i162 - 1][jxi468 - 1]);
dc0m[2 * dc0_index] = dc0_real;
dc0m[2 * dc0_index + 1] = dc0_imag;
dc0_index++;
@@ -926,8 +931,8 @@ void ScattererConfiguration::write_legacy(string file_name) {
if (i162 == 1) ici = ici + ies;
for (int i157 = 0; i157 < ici; i157++) {
double dc0_real, dc0_img;
- dc0_real = dc0_matrix[i157][i162 - 1][jxi468 - 1].real();
- dc0_img = dc0_matrix[i157][i162 - 1][jxi468 - 1].imag();
+ dc0_real = real(dc0_matrix[i157][i162 - 1][jxi468 - 1]);
+ dc0_img = imag(dc0_matrix[i157][i162 - 1][jxi468 - 1]);
// The FORTRAN code writes the complex numbers as a 16-byte long binary stream.
// Here we assume that the 16 bytes are equally split in 8 bytes to represent the
// real part and 8 bytes to represent the imaginary one.
@@ -1096,7 +1101,8 @@ void ScattererConfiguration::write_formatted(string file_name) {
if (i473 == 1) ici += ies;
fprintf(output, " SPHERE N.%4d\n", i473);
for (int ic472 = 0; ic472 < ici; ic472++) {
- double dc0_real = dc0_matrix[ic472][i473 - 1][0].real(), dc0_img = dc0_matrix[ic472][i473 - 1][0].imag();
+ double dc0_real = real(dc0_matrix[ic472][i473 - 1][0]);
+ double dc0_img = imag(dc0_matrix[ic472][i473 - 1][0]);
fprintf(output, "%5d %12.4lE%12.4lE\n", (ic472 + 1), dc0_real, dc0_img);
}
}
@@ -1110,8 +1116,8 @@ void ScattererConfiguration::write_formatted(string file_name) {
for (int ic477 = 1; ic477 <= ici; ic477++) {
fprintf(output, " NONTRANSITION LAYER N.%2d, SCALE = %3s\n", ic477, reference_variable_name.c_str());
for (int jxi476 = 0; jxi476 < number_of_scales; jxi476++) {
- double dc0_real = dc0_matrix[ic477 - 1][i478 - 1][jxi476].real();
- double dc0_img = dc0_matrix[ic477 - 1][i478 - 1][jxi476].imag();
+ double dc0_real = real(dc0_matrix[ic477 - 1][i478 - 1][jxi476]);
+ double dc0_img = imag(dc0_matrix[ic477 - 1][i478 - 1][jxi476]);
fprintf(output, "%5d %12.4lE%12.4lE\n", (jxi476 + 1), dc0_real, dc0_img);
}
}
diff --git a/src/libnptm/Makefile b/src/libnptm/Makefile
index b2fbc322c6b1a227ba0e29e80cd2d8e8c2c688bd..43b886ffe2745399f4e1389471c4050f1d0c4c48 100644
--- a/src/libnptm/Makefile
+++ b/src/libnptm/Makefile
@@ -19,11 +19,11 @@ endif
include ../make.inc
-CXX_NPTM_OBJS=$(OBJDIR)/Commons.o $(OBJDIR)/Configuration.o $(OBJDIR)/file_io.o $(OBJDIR)/Parsers.o $(OBJDIR)/sph_subs.o $(OBJDIR)/clu_subs.o $(OBJDIR)/tfrfme.o $(OBJDIR)/tra_subs.o $(OBJDIR)/TransitionMatrix.o
+CXX_NPTM_OBJS=$(OBJDIR)/Commons.o $(OBJDIR)/Configuration.o $(OBJDIR)/file_io.o $(OBJDIR)/Parsers.o $(OBJDIR)/sph_subs.o $(OBJDIR)/clu_subs.o $(OBJDIR)/tfrfme.o $(OBJDIR)/tra_subs.o $(OBJDIR)/TransitionMatrix.o $(OBJDIR)/lapack_calls.o $(OBJDIR)/algebraic.o $(OBJDIR)/types.o
-CXX_NPTM_DYNOBJS=$(DYNOBJDIR)/Commons.o $(DYNOBJDIR)/Configuration.o $(DYNOBJDIR)/file_io.o $(DYNOBJDIR)/Parsers.o $(DYNOBJDIR)/sph_subs.o $(DYNOBJDIR)/clu_subs.o $(DYNOBJDIR)/tfrfme.o $(DYNOBJDIR)/tra_subs.o $(DYNOBJDIR)/TransitionMatrix.o
+CXX_NPTM_DYNOBJS=$(DYNOBJDIR)/Commons.o $(DYNOBJDIR)/Configuration.o $(DYNOBJDIR)/file_io.o $(DYNOBJDIR)/Parsers.o $(DYNOBJDIR)/sph_subs.o $(DYNOBJDIR)/clu_subs.o $(DYNOBJDIR)/tfrfme.o $(DYNOBJDIR)/tra_subs.o $(DYNOBJDIR)/TransitionMatrix.o $(DYNOBJDIR)/lapack_calls.o $(DYNOBJDIR)/algebraic.o $(DYNOBJDIR)/types.o
-CXX_NPTM_DEBUG=$(OBJDIR)/Commons.g* $(OBJDIR)/Configuration.g* $(OBJDIR)/file_io.g* $(OBJDIR)/Parsers.g* $(OBJDIR)/sph_subs.g* $(OBJDIR)/clu_subs.g* $(OBJDIR)/tfrfme.g* $(OBJDIR)/tra_subs.g* $(OBJDIR)/TransitionMatrix.g*
+CXX_NPTM_DEBUG=$(OBJDIR)/Commons.g* $(OBJDIR)/Configuration.g* $(OBJDIR)/file_io.g* $(OBJDIR)/Parsers.g* $(OBJDIR)/sph_subs.g* $(OBJDIR)/clu_subs.g* $(OBJDIR)/tfrfme.g* $(OBJDIR)/tra_subs.g* $(OBJDIR)/TransitionMatrix.g* $(OBJDIR)/lapack_calls.g* $(OBJDIR)/algebraic.g* $(OBJDIR)/types.g* $(DYNOBJDIR)/Commons.g* $(DYNOBJDIR)/Configuration.g* $(DYNOBJDIR)/file_io.g* $(DYNOBJDIR)/Parsers.g* $(DYNOBJDIR)/sph_subs.g* $(DYNOBJDIR)/clu_subs.g* $(DYNOBJDIR)/tfrfme.g* $(DYNOBJDIR)/tra_subs.g* $(DYNOBJDIR)/TransitionMatrix.g* $(DYNOBJDIR)/lapack_calls.g* $(DYNOBJDIR)/algebraic.g* $(DYNOBJDIR)/types.g*
all: $(BUILDDIR_NPTM)/libnptm.a $(BUILDDIR_NPTM)/libnptm.so
diff --git a/src/libnptm/TransitionMatrix.cpp b/src/libnptm/TransitionMatrix.cpp
index 30126ff17ad8cf8511ce395ce3fa0df3f69b80e4..fc6c92af6c77dca0aabce667de72f1e51ef15fec 100644
--- a/src/libnptm/TransitionMatrix.cpp
+++ b/src/libnptm/TransitionMatrix.cpp
@@ -1,13 +1,18 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file TransitionMatrix.cpp
*
* \brief Implementation of the Transition Matrix structure.
*/
-#include <complex>
#include <exception>
#include <fstream>
#include <hdf5.h>
#include <string>
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
+
#ifndef INCLUDE_ERRORS_H_
#include "../include/errors.h"
#endif
@@ -32,7 +37,7 @@ TransitionMatrix::~TransitionMatrix() {
}
TransitionMatrix::TransitionMatrix(
- int _is, int _lm, double _vk, double _exri, std::complex<double> *_elements,
+ int _is, int _lm, double _vk, double _exri, dcomplex *_elements,
double _radius
) {
is = _is;
@@ -53,8 +58,8 @@ TransitionMatrix::TransitionMatrix(
}
TransitionMatrix::TransitionMatrix(
- int _lm, double _vk, double _exri, std::complex<double> **_rmi,
- std::complex<double> **_rei, double _sphere_radius
+ int _lm, double _vk, double _exri, dcomplex **_rmi,
+ dcomplex **_rei, double _sphere_radius
) {
is = 1111;
shape = new int[2];
@@ -64,7 +69,7 @@ TransitionMatrix::TransitionMatrix(
vk = _vk;
exri = _exri;
sphere_radius = _sphere_radius;
- elements = new complex<double>[2 * l_max]();
+ elements = new dcomplex[2 * l_max]();
for (int ei = 0; ei < l_max; ei++) {
elements[2 * ei] = -1.0 / _rmi[ei][0];
elements[2 * ei + 1] = -1.0 / _rei[ei][0];
@@ -73,7 +78,7 @@ TransitionMatrix::TransitionMatrix(
TransitionMatrix::TransitionMatrix(
int _nlemt, int _lm, double _vk, double _exri,
- std::complex<double> **_am0m
+ dcomplex **_am0m
) {
is = 1;
shape = new int[2];
@@ -83,7 +88,7 @@ TransitionMatrix::TransitionMatrix(
vk = _vk;
exri = _exri;
sphere_radius = 0.0;
- elements = new complex<double>[_nlemt * _nlemt]();
+ elements = new dcomplex[_nlemt * _nlemt]();
for (int ei = 0; ei < _nlemt; ei++) {
for (int ej = 0; ej < _nlemt; ej++) elements[_nlemt * ei + ej] = _am0m[ei][ej];
}
@@ -113,7 +118,7 @@ TransitionMatrix* TransitionMatrix::from_hdf5(string file_name) {
double _vk;
double _exri;
// This vector will be passed to the new object. DO NOT DELETE HERE!
- complex<double> *_elements;
+ dcomplex *_elements;
double _radius = 0.0;
status = hdf_file->read("IS", "INT32", &_is);
status = hdf_file->read("L_MAX", "INT32", &_lm);
@@ -125,9 +130,9 @@ TransitionMatrix* TransitionMatrix::from_hdf5(string file_name) {
hid_t file_id = hdf_file->get_file_id();
hid_t dset_id = H5Dopen2(file_id, "ELEMENTS", H5P_DEFAULT);
status = H5Dread(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, file_vector);
- _elements = new complex<double>[num_elements]();
+ _elements = new dcomplex[num_elements]();
for (int ei = 0; ei < num_elements; ei++) {
- _elements[ei] = complex<double>(file_vector[2 * ei], file_vector[2 * ei + 1]);
+ _elements[ei] = file_vector[2 * ei] + file_vector[2 * ei + 1] * I;
}
status = H5Dclose(dset_id);
status = hdf_file->read("RADIUS", "FLOAT64", &_radius);
@@ -140,9 +145,9 @@ TransitionMatrix* TransitionMatrix::from_hdf5(string file_name) {
hid_t file_id = hdf_file->get_file_id();
hid_t dset_id = H5Dopen2(file_id, "ELEMENTS", H5P_DEFAULT);
status = H5Dread(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, file_vector);
- _elements = new complex<double>[num_elements]();
+ _elements = new dcomplex[num_elements]();
for (int ei = 0; ei < num_elements; ei++) {
- _elements[ei] = complex<double>(file_vector[2 * ei], file_vector[2 * ei + 1]);
+ _elements[ei] = file_vector[2 * ei] + file_vector[2 * ei + 1] * I;
}
status = H5Dclose(dset_id);
tm = new TransitionMatrix(_is, _lm, _vk, _exri, _elements, _radius);
@@ -166,7 +171,7 @@ TransitionMatrix* TransitionMatrix::from_legacy(string file_name) {
double _vk;
double _exri;
// This vector will be passed to the new object. DO NOT DELETE HERE!
- complex<double> *_elements;
+ dcomplex *_elements;
double _radius = 0.0;
ttms.read(reinterpret_cast<char *>(&_is), sizeof(int));
ttms.read(reinterpret_cast<char *>(&_lm), sizeof(int));
@@ -174,12 +179,12 @@ TransitionMatrix* TransitionMatrix::from_legacy(string file_name) {
ttms.read(reinterpret_cast<char *>(&_exri), sizeof(double));
if (_is == 1111) {
num_elements = _lm * 2;
- _elements = new complex<double>[num_elements]();
+ _elements = new dcomplex[num_elements]();
for (int ei = 0; ei < num_elements; ei++) {
double vreal, vimag;
ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));
ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
- _elements[ei] = complex<double>(vreal, vimag);
+ _elements[ei] = vreal + vimag * I;
}
double _radius;
ttms.read(reinterpret_cast<char *>(&_radius), sizeof(double));
@@ -187,12 +192,12 @@ TransitionMatrix* TransitionMatrix::from_legacy(string file_name) {
} else if (_is == 1) {
int nlemt = 2 * _lm * (_lm + 2);
num_elements = nlemt * nlemt;
- _elements = new complex<double>[num_elements]();
+ _elements = new dcomplex[num_elements]();
for (int ei = 0; ei < num_elements; ei++) {
double vreal, vimag;
ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));
ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
- _elements[ei] = complex<double>(vreal, vimag);
+ _elements[ei] = vreal + vimag * I;
}
tm = new TransitionMatrix(_is, _lm, _vk, _exri, _elements);
}
@@ -239,8 +244,8 @@ void TransitionMatrix::write_hdf5(string file_name) {
int num_elements = 2 * shape[0] * shape[1];
double *ptr_elements = new double[num_elements]();
for (int ei = 0; ei < num_elements / 2; ei++) {
- ptr_elements[2 * ei] = elements[ei].real();
- ptr_elements[2 * ei + 1] = elements[ei].imag();
+ ptr_elements[2 * ei] = real(elements[ei]);
+ ptr_elements[2 * ei + 1] = imag(elements[ei]);
}
rec_ptr_list.append(ptr_elements);
if (is == 1111) {
@@ -287,10 +292,10 @@ void TransitionMatrix::write_legacy(string file_name) {
if (ttms.is_open()) {
int num_elements = shape[0] * shape[1];
for (int ei = 0; ei < num_elements; ei++) {
- complex<double> element1 = elements[ei];
+ dcomplex element1 = elements[ei];
double vreal, vimag;
- vreal = element1.real();
- vimag = element1.imag();
+ vreal = real(element1);
+ vimag = imag(element1);
ttms.write(reinterpret_cast<char *>(&vreal), sizeof(double));
ttms.write(reinterpret_cast<char *>(&vimag), sizeof(double));
}
diff --git a/src/libnptm/algebraic.cpp b/src/libnptm/algebraic.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..48335792275ba7097b29f5b639fe828bbd1aa8d2
--- /dev/null
+++ b/src/libnptm/algebraic.cpp
@@ -0,0 +1,39 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
+/*! \file algebraic.cpp
+ *
+ * \brief Implementation of algebraic functions with different call-backs.
+ */
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
+
+#ifdef USE_LAPACK
+#ifdef USE_MKL
+#include <mkl_lapacke.h>
+#else
+#include <lapacke.h>
+#endif
+#ifndef INCLUDE_LAPACK_CALLS_H_
+#include "../include/lapack_calls.h"
+#endif
+#endif
+
+#ifndef INCLUDE_ALGEBRAIC_H_
+#include "../include/algebraic.h"
+#endif
+
+// >>> FALL-BACK FUNCTIONS DECLARATION <<< //
+extern void lucin(dcomplex **mat, np_int max_size, np_int size, int &ier);
+// >>> END OF FALL-BACK FUNCTIONS <<< //
+
+using namespace std;
+
+void invert_matrix(dcomplex **mat, np_int size, int &ier, np_int max_size) {
+ ier = 0;
+#ifdef USE_LAPACK
+ zinvert(mat, size, ier);
+#else
+ lucin(mat, max_size, size, ier);
+#endif
+}
diff --git a/src/libnptm/clu_subs.cpp b/src/libnptm/clu_subs.cpp
index 22287bb7c734fcfc067c4e185b4b4f7932c03918..72cac9914a8d1d61a6f67e93783871476f03ea21 100644
--- a/src/libnptm/clu_subs.cpp
+++ b/src/libnptm/clu_subs.cpp
@@ -1,8 +1,13 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file clu_subs.cpp
*
* \brief C++ implementation of CLUSTER subroutines.
*/
-#include <complex>
+
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
#ifndef INCLUDE_COMMONS_H_
#include "../include/Commons.h"
@@ -19,22 +24,22 @@
using namespace std;
void apc(
- double ****zpv, int le, std::complex<double> **am0m, std::complex<double> **w,
- double sqk, double **gapr, std::complex<double> **gapp
+ double ****zpv, int le, dcomplex **am0m, dcomplex **w,
+ double sqk, double **gapr, dcomplex **gapp
) {
- complex<double> **ac, **gap;
- const complex<double> cc0(0.0, 0.0);
- const complex<double> uim(0.0, 1.0);
- complex<double> uimmp, summ, sume, suem, suee, summp, sumep;
- complex<double> suemp, sueep;
+ dcomplex **ac, **gap;
+ const dcomplex cc0 = 0.0 + 0.0 * I;
+ const dcomplex uim = 0.0 + 1.0 * I;
+ dcomplex uimmp, summ, sume, suem, suee, summp, sumep;
+ dcomplex suemp, sueep;
double cof = 1.0 / sqk;
double cimu = cof / sqrt(2.0);
int nlem = le * (le + 2);
const int nlemt = nlem + nlem;
- ac = new complex<double>*[nlemt];
- gap = new complex<double>*[3];
- for (int ai = 0; ai < nlemt; ai++) ac[ai] = new complex<double>[2]();
- for (int gi = 0; gi < 3; gi++) gap[gi] = new complex<double>[2]();
+ ac = new dcomplex*[nlemt];
+ gap = new dcomplex*[3];
+ for (int ai = 0; ai < nlemt; ai++) ac[ai] = new dcomplex[2]();
+ for (int gi = 0; gi < 3; gi++) gap[gi] = new dcomplex[2]();
for (int j45 = 1; j45 <= nlemt; j45++) {
int j = j45 - 1;
ac[j][0] = cc0;
@@ -119,9 +124,9 @@ void apc(
sume = gap[0][ipo95 - 1] * cimu;
suee = gap[1][ipo95 - 1] * cof;
suem = gap[2][ipo95 - 1] * cimu;
- gapr[0][ipo95 - 1] = (sume - suem).real();
- gapr[1][ipo95 - 1] = ((sume + suem) * uim).real();
- gapr[2][ipo95 - 1] = suee.real();
+ gapr[0][ipo95 - 1] = real(sume - suem);
+ gapr[1][ipo95 - 1] = real((sume + suem) * uim);
+ gapr[2][ipo95 - 1] = real(suee);
sumep = gapp[0][ipo95 - 1] * cimu;
sueep = gapp[1][ipo95 - 1] * cof;
suemp = gapp[2][ipo95 - 1] * cimu;
@@ -137,24 +142,24 @@ void apc(
}
void apcra(
- double ****zpv, const int le, std::complex<double> **am0m, int inpol, double sqk,
- double **gaprm, std::complex<double> **gappm
+ double ****zpv, const int le, dcomplex **am0m, int inpol, double sqk,
+ double **gaprm, dcomplex **gappm
) {
- const complex<double> cc0(0.0, 0.0);
- const complex<double> uim(0.0, 1.0);
- complex<double> uimtl, uimtls, ca11, ca12, ca21, ca22;
- complex<double> a11, a12, a21, a22, sum1, sum2, fc;
+ const dcomplex cc0 = 0.0 + 0.0 * I;
+ const dcomplex uim = 0.0 + 1.0 * I;
+ dcomplex uimtl, uimtls, ca11, ca12, ca21, ca22;
+ dcomplex a11, a12, a21, a22, sum1, sum2, fc;
double ****svw = new double***[le];
- complex<double> ****svs = new complex<double>***[le];
+ dcomplex ****svs = new dcomplex***[le];
for (int i = 0; i < le; i++) {
svw[i] = new double**[3];
- svs[i] = new complex<double>**[3];
+ svs[i] = new dcomplex**[3];
for (int j = 0; j < 3; j++) {
svw[i][j] = new double*[2];
- svs[i][j] = new complex<double>*[2];
+ svs[i][j] = new dcomplex*[2];
for (int k = 0; k < 2; k++) {
svw[i][j][k] = new double[2]();
- svs[i][j][k] = new complex<double>[2]();
+ svs[i][j][k] = new dcomplex[2]();
}
}
}
@@ -206,7 +211,7 @@ void apcra(
int lmp = l58 + lmpml;
int impmm = lmp * (lmp + 1);
uimtl = uim * (1.0 * lmpml);
- if (lmpml == 0) uimtl = complex<double>(1.0, 0.0);
+ if (lmpml == 0) uimtl = 1.0 + 0.0 * I;
for (int im54 = 1; im54 <= ltpo; im54++) {
int m = im54 - lpo;
int i = imm + m;
@@ -229,7 +234,7 @@ void apcra(
int lsmp = ls + lsmpml;
int ismpmm = lsmp * (lsmp + 1);
uimtls = -uim * (1.0 * lsmpml);
- if (lsmpml == 0) uimtls = complex<double>(1.0, 0.0);
+ if (lsmpml == 0) uimtls = 1.0 + 0.0 * I;
for (int ims = 1; ims <= lstpo; ims++) {
int ms = ims - lspo;
int msmp = ms - mu;
@@ -328,14 +333,14 @@ void apcra(
if (inpol == 0) {
sum1 *= cofs;
sum2 *= cofs;
- gaprm[2][0] = sum1.real();
- gaprm[2][1] = sum1.real();
+ gaprm[2][0] = real(sum1);
+ gaprm[2][1] = real(sum1);
gappm[2][0] = sum2 * uim;
gappm[2][1] = -gappm[2][0];
} else { // label 72
cofs *= 2.0;
- gaprm[2][0] = sum1.real() * cofs;
- gaprm[2][1] = sum2.real() * cofs;
+ gaprm[2][0] = real(sum1) * cofs;
+ gaprm[2][1] = real(sum2) * cofs;
gappm[2][0] = cc0;
gappm[2][1] = cc0;
}
@@ -357,15 +362,12 @@ void apcra(
delete[] svs;
}
-complex<double> cdtp(
- std::complex<double> z, std::complex<double> **am, int i, int jf,
- int k, int nj
-) {
+dcomplex cdtp(dcomplex z, dcomplex **am, int i, int jf, int k, int nj) {
/* NOTE: the original FORTRAN code treats the AM matrix as a
* vector. This is not directly allowed in C++ and it requires
* accounting for the different dimensions.
*/
- complex<double> result = z;
+ dcomplex result = z;
if (nj > 0) {
int jl = jf + nj - 1;
for (int j = jf; j <= jl; j++) {
@@ -409,9 +411,9 @@ double cgev(int ipamo, int mu, int l, int m) {
return result;
}
-void cms(std::complex<double> **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6) {
- complex<double> dm, de, cgh, cgk;
- const complex<double> cc0(0.0, 0.0);
+void cms(dcomplex **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6) {
+ dcomplex dm, de, cgh, cgk;
+ const dcomplex cc0 = 0.0 + 0.0 * I;
int ndi = c4->nsph * c4->nlim;
int nbl = 0;
int nsphmo = c4->nsph - 1;
@@ -492,19 +494,19 @@ void cms(std::complex<double> **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6) {
}
void crsm1(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6) {
- complex<double> ***svf, ***svw, **svs;
- const complex<double> cc0(0.0, 0.0);
- complex<double> cam(0.0, 0.0);
+ dcomplex ***svf, ***svw, **svs;
+ const dcomplex cc0 = 0.0 + 0.0 * I;
+ dcomplex cam = cc0;
const int le4po = 4 * c4->le + 1;
- svf = new complex<double>**[le4po];
- svw = new complex<double>**[le4po];
- svs = new complex<double>*[le4po];
+ svf = new dcomplex**[le4po];
+ svw = new dcomplex**[le4po];
+ svs = new dcomplex*[le4po];
for (int si = 0; si < le4po; si++) {
- svf[si] = new complex<double>*[le4po];
- svw[si] = new complex<double>*[4];
- svs[si] = new complex<double>[4]();
- for (int sj = 0; sj < le4po; sj++) svf[si][sj] = new complex<double>[4]();
- for (int sj = 0; sj < 4; sj++) svw[si][sj] = new complex<double>[4]();
+ svf[si] = new dcomplex*[le4po];
+ svw[si] = new dcomplex*[4];
+ svs[si] = new dcomplex[4]();
+ for (int sj = 0; sj < le4po; sj++) svf[si][sj] = new dcomplex[4]();
+ for (int sj = 0; sj < 4; sj++) svw[si][sj] = new dcomplex[4]();
}
double exdc = exri * exri;
double ccs = 1.0 / (vk * vk);
@@ -611,22 +613,22 @@ void crsm1(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6) {
delete[] svs;
}
-complex<double> ghit(
- int ihi, int ipamo, int nbl, int l1, int m1, int l2, int m2, C1 *c1,
- C1_AddOns *c1ao, C4 *c4, C6 *c6
+dcomplex ghit(
+ int ihi, int ipamo, int nbl, int l1, int m1, int l2, int m2, C1 *c1,
+ C1_AddOns *c1ao, C4 *c4, C6 *c6
) {
/* NBL identifies transfer vector going from N2 to N1;
* IHI=0 for Hankel, IHI=1 for Bessel, IHI=2 for Bessel from origin;
* depending on IHI, IPAM=0 gives H or I, IPAM= 1 gives K or L. */
- const complex<double> cc0(0.0, 0.0);
- const complex<double> uim(0.0, 1.0);
- complex<double> csum(0.0, 0.0), cfun(0.0, 0.0);
- complex<double> result = cc0;
+ const dcomplex cc0 = 0.0 + 0.0 * I;
+ const dcomplex uim = 0.0 + 1.0 * I;
+ dcomplex csum = cc0, cfun = cc0;
+ dcomplex result = cc0;
if (ihi == 2) {
if (c1->rxx[nbl - 1] == 0.0 && c1->ryy[nbl - 1] == 0.0 && c1->rzz[nbl - 1] == 0.0) {
if (ipamo == 0) {
- if (l1 == l2 && m1 == m2) result = complex<double>(1.0, 0.0);
+ if (l1 == l2 && m1 == m2) result = 1.0 + 0.0 * I;
}
return result;
}
@@ -823,7 +825,7 @@ complex<double> ghit(
}
void hjv(
- double exri, double vk, int &jer, int &lcalc, std::complex<double> &arg,
+ double exri, double vk, int &jer, int &lcalc, dcomplex &arg,
C1 *c1, C1_AddOns *c1ao, C4 *c4
) {
int nsphmo = c4->nsph - 1;
@@ -844,7 +846,7 @@ void hjv(
double rz = c1->rzz[nf40 - 1] - c1->rzz[ns40 - 1];
double rr = sqrt(rx * rx + ry * ry + rz * rz);
double rarg = rr * vk * exri;
- arg = complex<double>(rarg, 0.0);
+ arg = rarg + 0.0 * I;
rbf(lit, rarg, lcalc, rfj);
if (lcalc < lit) {
jer = 1;
@@ -862,7 +864,7 @@ void hjv(
for (int lpo38 = 1; lpo38 <= c4->litpo; lpo38++) {
double rpart = rfj[lpo38 - 1];
double ipart = rfn[lpo38 - 1];
- c1ao->vh[lpo38 + ivhb - 1] = complex<double>(rpart, ipart);
+ c1ao->vh[lpo38 + ivhb - 1] = rpart + ipart * I;
}
ivhb += c4->litpo;
} // ns40 loop
@@ -892,23 +894,23 @@ void hjv(
delete[] rfn;
}
-void lucin(std::complex<double> **am, const int nddmst, int n, int &ier) {
+void lucin(dcomplex **am, const np_int nddmst, np_int n, int &ier) {
/* NDDMST FIRST DIMENSION OF AM AS DECLARED IN DIMENSION
* STATEMENT.
* N NUMBER OF ROWS IN AM.
* IER IS REPLACED BY 1 FOR SINGULARITY.
*/
double *v = new double[nddmst];
- complex<double> ctemp, cfun;
- complex<double> cc0 = complex<double>(0.0, 0.0);
+ dcomplex ctemp, cfun;
+ const dcomplex cc0 = 0.0 + 0.0 * I;
ier = 0;
int nminus = n - 1;
- for (int i = 1; i <= n; i++) {
+ for (int64_t i = 1; i <= n; i++) {
double sum = 0.0;
- for (int j = 1; j <= n; j++) {
+ for (int64_t j = 1; j <= n; j++) {
sum += (
- am[i - 1][j - 1].real() * am[i - 1][j - 1].real()
- + am[i - 1][j - 1].imag() * am[i - 1][j - 1].imag()
+ real(am[i - 1][j - 1]) * real(am[i - 1][j - 1])
+ + imag(am[i - 1][j - 1]) * imag(am[i - 1][j - 1])
);
} // j1319 loop
v[i - 1] = 1.0 / sum;
@@ -917,24 +919,23 @@ void lucin(std::complex<double> **am, const int nddmst, int n, int &ier) {
// REPLACE L(I,I) BY 1/L(I,I), READY FOR SECTION 4.
// (ROW INTERCHANGES TAKE PLACE, AND THE INDICES OF THE PIVOTAL ROWS
// ARE PLACED IN V.)
- /* >>> THERE APPEARS TO BE A BUG IN THE FOLLOWING LOOP <<< */
- for (int k = 1; k <= n; k++) {
- int kplus = k + 1;
- int kminus = k - 1;
- int l = k;
+ for (int64_t k = 1; k <= n; k++) {
+ int64_t kplus = k + 1;
+ int64_t kminus = k - 1;
+ int64_t l = k;
double psqmax = 0.0;
- for (int i = k; i <= n; i++) {
+ for (int64_t i = k; i <= n; i++) {
cfun = cdtp(-am[i - 1][k - 1], am, i, 1, k, kminus);
ctemp = -cfun;
am[i - 1][k - 1] = ctemp;
- double psq = v[i - 1] * (ctemp.real() * ctemp.real() + ctemp.imag() * ctemp.imag());
+ double psq = v[i - 1] * (real(ctemp) * real(ctemp) + imag(ctemp) * imag(ctemp));
if (psq > psqmax) {
psqmax = psq;
l = i;
}
} // i2029 loop
if (l != k) {
- for (int j = 1; j <= n; j++) {
+ for (int64_t j = 1; j <= n; j++) {
ctemp = am[k - 1][j - 1];
am[k - 1][j - 1] = am[l - 1][j - 1];
am[l - 1][j - 1] = ctemp;
@@ -951,7 +952,7 @@ void lucin(std::complex<double> **am, const int nddmst, int n, int &ier) {
ctemp = 1.0 / am[k - 1][k - 1];
am[k - 1][k - 1] = ctemp;
if (kplus <= n) {
- for (int j = kplus; j <= n; j++) {
+ for (int64_t j = kplus; j <= n; j++) {
cfun = cdtp(-am[k - 1][j - 1], am, k, 1, j, kminus);
am[k - 1][j - 1] = -ctemp * cfun;
} // j2059 loop
@@ -959,9 +960,9 @@ void lucin(std::complex<double> **am, const int nddmst, int n, int &ier) {
} // k2019 loop
// 4. REPLACE AM BY ITS INVERSE AMINV.
// 4.1 REPLACE L AND U BY THEIR INVERSE LINV AND UINV.
- for (int k = 1; k <= nminus; k++) {
- int kplus = k + 1;
- for (int i = kplus; i <= n; i++) {
+ for (int64_t k = 1; k <= nminus; k++) {
+ int64_t kplus = k + 1;
+ for (int64_t i = kplus; i <= n; i++) {
cfun = cdtp(cc0, am, i, k, k, i - k);
am[i - 1][k - 1] = -am[i - 1][i - 1] * cfun;
cfun = cdtp(am[k - 1][i - 1], am, k, kplus, i, i - k - 1);
@@ -969,8 +970,8 @@ void lucin(std::complex<double> **am, const int nddmst, int n, int &ier) {
} // i4119 loop
} // k4109 loop
// 4.2 FORM AMINV=UINV*LINV.
- for (int k = 1; k <= n; k++) {
- for (int i = 1; i <= n; i++) {
+ for (int64_t k = 1; k <= n; k++) {
+ for (int64_t i = 1; i <= n; i++) {
if (i < k) {
cfun = cdtp(cc0, am, i, k, k, n - k + 1);
am[i - 1][k -1] = cfun;
@@ -983,11 +984,11 @@ void lucin(std::complex<double> **am, const int nddmst, int n, int &ier) {
} // k4209 loop
// 4.3 INTERCHANGE COLUMNS OF AMINV AS SPECIFIED BY V, BUT IN REVERSE
// ORDER.
- for (int l = 1; l <= n; l++) {
- int k = n - l + 1;
- int kcol = (int)(v[k - 1]);
+ for (int64_t l = 1; l <= n; l++) {
+ int64_t k = n - l + 1;
+ int64_t kcol = (int64_t)(v[k - 1]);
if (kcol != k) {
- for (int i = 1; i <= n; i++) {
+ for (int64_t i = 1; i <= n; i++) {
ctemp = am[i - 1][k - 1];
am[i - 1][k - 1] = am[i - 1][kcol - 1];
am[i - 1][kcol - 1] = ctemp;
@@ -997,15 +998,15 @@ void lucin(std::complex<double> **am, const int nddmst, int n, int &ier) {
delete[] v;
}
-void mextc(double vk, double exri, std::complex<double> **fsac, double **cextlr, double **cext) {
- double fa11r = fsac[0][0].real();
- double fa11i = fsac[0][0].imag();
- double fa21r = fsac[1][0].real();
- double fa21i = fsac[1][0].imag();
- double fa12r = fsac[0][1].real();
- double fa12i = fsac[0][1].imag();
- double fa22r = fsac[1][1].real();
- double fa22i = fsac[1][1].imag();
+void mextc(double vk, double exri, dcomplex **fsac, double **cextlr, double **cext) {
+ double fa11r = real(fsac[0][0]);
+ double fa11i = imag(fsac[0][0]);
+ double fa21r = real(fsac[1][0]);
+ double fa21i = imag(fsac[1][0]);
+ double fa12r = real(fsac[0][1]);
+ double fa12i = imag(fsac[0][1]);
+ double fa22r = real(fsac[1][1]);
+ double fa22i = imag(fsac[1][1]);
cextlr[0][0] = fa11i * 2.0;
cextlr[0][1] = 0.0;
cextlr[0][2] = -fa12i;
@@ -1048,12 +1049,12 @@ void mextc(double vk, double exri, std::complex<double> **fsac, double **cextlr,
}
void pcros(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C4 *c4) {
- const complex<double> cc0(0.0, 0.0);
- complex<double> sump, sum1, sum2, sum3, sum4, am, amp, cc, csam;
+ const dcomplex cc0 = 0.0 + 0.0 * I;
+ dcomplex sump, sum1, sum2, sum3, sum4, am, amp, cc, csam;
const double exdc = exri * exri;
double ccs = 1.0 / (vk * vk);
double cccs = ccs / exdc;
- csam = -(ccs / (exri * vk)) * complex<double>(0.0, 0.5);
+ csam = -(ccs / (exri * vk)) * 0.5 * I;
const double pi4sq = 64.0 * acos(0.0) * acos(0.0);
double cfsq = 4.0 / (pi4sq *ccs * ccs);
const int nlemt = c4->nlem + c4->nlem;
@@ -1077,7 +1078,7 @@ void pcros(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C4 *c4) {
am += (c1ao->am0m[i][j] * c1->w[j][ipo18 - 1]);
amp += (c1ao->am0m[i][j] * c1->w[j][jpo - 1]);
} // j10 loop
- sum += (dconjg(am) * am).real();
+ sum += real(dconjg(am) * am);
sump += (dconjg(amp) * am);
sum1 += (dconjg(c1->w[i][ipo18 - 1]) * am);
sum2 += (dconjg(c1->w[i][jpo - 1]) * am);
@@ -1086,7 +1087,7 @@ void pcros(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C4 *c4) {
} // i12 loop
c1ao->scsc[ipo18 - 1] = cccs * sum;
c1ao->scscp[ipo18 - 1] = cccs * sump;
- c1ao->ecsc[ipo18 - 1] = -cccs * sum1.real();
+ c1ao->ecsc[ipo18 - 1] = -cccs * real(sum1);
c1ao->ecscp[ipo18 - 1] = -cccs * sum2;
c1ao->fsac[ipo18 - 1][ipo18 - 1] = csam * sum1;
c1ao->fsac[jpo - 1][ipo18 - 1] = csam * sum2;
@@ -1107,14 +1108,14 @@ void pcros(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C4 *c4) {
}
void pcrsm0(double vk, double exri, int inpol, C1 *c1, C1_AddOns *c1ao, C4 *c4) {
- const complex<double> cc0(0.0, 0.0);
- const complex<double> uim(0.0, 1.0);
- complex<double> sum1, sum2, sum3, sum4, sumpd;
- complex<double> sums1, sums2, sums3, sums4, csam;
+ const dcomplex cc0 = 0.0 + 0.0 * I;
+ const dcomplex uim = 0.0 + 1.0 * I;
+ dcomplex sum1, sum2, sum3, sum4, sumpd;
+ dcomplex sums1, sums2, sums3, sums4, csam;
double exdc = exri * exri;
double ccs = 4.0 * acos(0.0) / (vk * vk);
double cccs = ccs / exdc;
- csam = -(ccs / (exri * vk)) * complex<double>(0.0, 0.5);
+ csam = -(ccs / (exri * vk)) * 0.5 * I;
sum2 = cc0;
sum3 = cc0;
for (int i14 = 1; i14 <= c4->nlem; i14++) { // GPU portable?
@@ -1128,10 +1129,10 @@ void pcrsm0(double vk, double exri, int inpol, C1 *c1, C1_AddOns *c1ao, C4 *c4)
for (int i16 = 1; i16 <= nlemt; i16++) {
for (int j16 = 1; j16 <= c4->nlem; j16++) {
int je = j16 + c4->nlem;
- double rvalue = (
- dconjg(c1ao->am0m[i16 - 1][j16 - 1]) * c1ao->am0m[i16 - 1][j16 - 1]
- + dconjg(c1ao->am0m[i16 - 1][je - 1]) * c1ao->am0m[i16 - 1][je - 1]
- ).real();
+ double rvalue = real(
+ dconjg(c1ao->am0m[i16 - 1][j16 - 1]) * c1ao->am0m[i16 - 1][j16 - 1]
+ + dconjg(c1ao->am0m[i16 - 1][je - 1]) * c1ao->am0m[i16 - 1][je - 1]
+ );
sumpi += rvalue;
sumpd += (
dconjg(c1ao->am0m[i16 - 1][j16 - 1]) * c1ao->am0m[i16 - 1][je - 1]
@@ -1158,16 +1159,16 @@ void pcrsm0(double vk, double exri, int inpol, C1 *c1, C1_AddOns *c1ao, C4 *c4)
sums4 = cc0;
}
// label 20
- c1ao->ecscm[0] = -cccs * sum2.real();
- c1ao->ecscm[1] = -cccs * sum1.real();
+ c1ao->ecscm[0] = -cccs * real(sum2);
+ c1ao->ecscm[1] = -cccs * real(sum1);
c1ao->ecscpm[0] = -cccs * sum4;
c1ao->ecscpm[1] = -cccs * sum3;
c1ao->fsacm[0][0] = csam * sum2;
c1ao->fsacm[1][0] = csam * sum4;
c1ao->fsacm[1][1] = csam * sum1;
c1ao->fsacm[0][1] = csam * sum3;
- c1ao->scscm[0] = cccs * sums1.real();
- c1ao->scscm[1] = cccs * sums2.real();
+ c1ao->scscm[0] = cccs * real(sums1);
+ c1ao->scscm[1] = cccs * real(sums2);
c1ao->scscpm[0] = cccs * sums3;
c1ao->scscpm[1] = cccs * sums4;
}
@@ -1538,23 +1539,23 @@ void r3jmr(int j1, int j2, int j3, int m1, C6 *c6) {
}
void raba(
- int le, std::complex<double> **am0m, std::complex<double> **w, double **tqce,
- std::complex<double> **tqcpe, double **tqcs, std::complex<double> **tqcps
+ int le, dcomplex **am0m, dcomplex **w, double **tqce,
+ dcomplex **tqcpe, double **tqcs, dcomplex **tqcps
) {
- complex<double> **a, **ctqce, **ctqcs;
- complex<double> acw, acwp, aca, acap, c1, c2, c3;
- const complex<double> cc0(0.0, 0.0);
- const complex<double> uim(0.0, 1.0);
+ dcomplex **a, **ctqce, **ctqcs;
+ dcomplex acw, acwp, aca, acap, c1, c2, c3;
+ const dcomplex cc0 = 0.0 + 0.0 * I;
+ const dcomplex uim = 0.0 + 1.0 * I;
const double sq2i = 1.0 / sqrt(2.0);
int nlem = le * (le + 2);
const int nlemt = nlem + nlem;
- a = new complex<double>*[nlemt];
- ctqce = new complex<double>*[2];
- ctqcs = new complex<double>*[2];
- for (int ai = 0; ai < nlemt; ai++) a[ai] = new complex<double>[2]();
+ a = new dcomplex*[nlemt];
+ ctqce = new dcomplex*[2];
+ ctqcs = new dcomplex*[2];
+ for (int ai = 0; ai < nlemt; ai++) a[ai] = new dcomplex[2]();
for (int ci = 0; ci < 2; ci++) {
- ctqce[ci] = new complex<double>[3]();
- ctqcs[ci] = new complex<double>[3]();
+ ctqce[ci] = new dcomplex[3]();
+ ctqcs[ci] = new dcomplex[3]();
}
for (int i20 = 1; i20 <= nlemt; i20++) {
int i = i20 - 1;
@@ -1638,18 +1639,18 @@ void raba(
} // ipo70 loop
for (int ipo78 = 1; ipo78 <= 2; ipo78++) {
int ipo = ipo78 - 1;
- tqce[ipo][0] = (ctqce[ipo][0] - ctqce[ipo][2]).real() * sq2i;
- tqce[ipo][1] = ((ctqce[ipo][0] + ctqce[ipo][2]) * uim).real() * sq2i;
- tqce[ipo][2] = ctqce[ipo][1].real();
+ tqce[ipo][0] = real(ctqce[ipo][0] - ctqce[ipo][2]) * sq2i;
+ tqce[ipo][1] = real((ctqce[ipo][0] + ctqce[ipo][2]) * uim) * sq2i;
+ tqce[ipo][2] = real(ctqce[ipo][1]);
c1 = tqcpe[ipo][0];
c2 = tqcpe[ipo][1];
c3 = tqcpe[ipo][2];
tqcpe[ipo][0] = (c1 - c3) * sq2i;
tqcpe[ipo][1] = (c1 + c3) * (uim * sq2i);
tqcpe[ipo][2] = c2;
- tqcs[ipo][0] = -sq2i * (ctqcs[ipo][0] - ctqcs[ipo][2]).real();
- tqcs[ipo][1] = -sq2i * ((ctqcs[ipo][0] + ctqcs[ipo][2]) * uim).real();
- tqcs[ipo][2] = -1.0 * ctqcs[ipo][1].real();
+ tqcs[ipo][0] = -sq2i * real(ctqcs[ipo][0] - ctqcs[ipo][2]);
+ tqcs[ipo][1] = -sq2i * real((ctqcs[ipo][0] + ctqcs[ipo][2]) * uim);
+ tqcs[ipo][2] = -1.0 * real(ctqcs[ipo][1]);
c1 = tqcps[ipo][0];
c2 = tqcps[ipo][1];
c3 = tqcps[ipo][2];
@@ -1685,12 +1686,12 @@ void rftr(
}
void scr0(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 * c4) {
- const complex<double> cc0(0.0, 0.0);
+ const dcomplex cc0 = 0.0 + 0.0 * I;
double exdc = exri * exri;
double ccs = 4.0 * acos(0.0) / (vk * vk);
double cccs = ccs / exdc;
- complex<double> sum21, rm, re, csam;
- csam = -(ccs / (exri * vk)) * complex<double>(0.0, 0.5);
+ dcomplex sum21, rm, re, csam;
+ csam = -(ccs / (exri * vk)) * 0.5 * I;
//double scs = 0.0, ecs = 0.0, acs = 0.0;
c3->scs = 0.0;
c3->ecs = 0.0;
@@ -1705,13 +1706,13 @@ void scr0(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 * c4) {
double fl = 1.0 * (l10 + l10 + 1);
rm = 1.0 / c1->rmi[l10 - 1][i14 - 1];
re = 1.0 / c1->rei[l10 - 1][i14 - 1];
- double rvalue = (dconjg(rm) * rm + dconjg(re) * re).real() * fl;
+ double rvalue = real(dconjg(rm) * rm + dconjg(re) * re) * fl;
sums += rvalue;
sum21 += ((rm + re) * fl);
} // l10 loop
sum21 *= -1.0;
double scasec = cccs * sums;
- double extsec = -cccs * sum21.real();
+ double extsec = -cccs * real(sum21);
double abssec = extsec - scasec;
c1->sscs[i14 - 1] = scasec;
c1->sexs[i14 - 1] = extsec;
@@ -1734,11 +1735,11 @@ void scr2(
double vk, double vkarg, double exri, double *duk, C1 *c1, C1_AddOns *c1ao,
C3 *c3, C4 *c4
) {
- const complex<double> cc0(0.0, 0.0);
- const complex<double> uim(0.0, 1.0);
- complex<double> s11, s21, s12, s22, rm, re, csam, cph, phas, cc;
+ const dcomplex cc0 = 0.0 + 0.0 * I;
+ const dcomplex uim = 0.0 + 1.0 * I;
+ dcomplex s11, s21, s12, s22, rm, re, csam, cph, phas, cc;
double ccs = 1.0 / (vk * vk);
- csam = -(ccs / (exri * vk)) * complex<double>(0.0, 0.5);
+ csam = -(ccs / (exri * vk)) * 0.5 * I;
const double pi4sq = 64.0 * acos(0.0) * acos(0.0);
double cfsq = 4.0 / (pi4sq * ccs * ccs);
cph = uim * exri * vkarg;
@@ -1776,7 +1777,7 @@ void scr2(
c1->sas[i][1][1] = s22 * csam;
}
// label 12
- phas = exp(cph * (duk[0] * c1->rxx[i] + duk[1] * c1->ryy[i] + duk[2] * c1->rzz[i]));
+ phas = cexp(cph * (duk[0] * c1->rxx[i] + duk[1] * c1->ryy[i] + duk[2] * c1->rzz[i]));
c3->tsas[0][0] += (c1->sas[iogi - 1][0][0] * phas);
c3->tsas[1][0] += (c1->sas[iogi - 1][1][0] * phas);
c3->tsas[0][1] += (c1->sas[iogi - 1][0][1] * phas);
@@ -1814,7 +1815,7 @@ void scr2(
}
void str(double **rcf, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 *c4, C6 *c6) {
- complex<double> *ylm;
+ dcomplex *ylm;
const double pi = acos(-1.0);
c3->gcs = 0.0;
double gcss = 0.0;
@@ -1831,7 +1832,7 @@ void str(double **rcf, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 *c4, C6 *c6) {
c3->gcs += gcss;
} // i18 loop
int ylm_size = (c4->litpos > c4->lmtpos) ? c4->litpos : c4->lmtpos;
- ylm = new complex<double>[ylm_size]();
+ ylm = new dcomplex[ylm_size]();
int i = 0;
for (int l1po28 = 1; l1po28 <= c4->lmpo; l1po28++) {
int l1 = l1po28 - 1;
@@ -1901,9 +1902,9 @@ void tqr(
tsk = u[0] * tqsv[0] + u[1] * tqsv[1] + u[2] * tqsv[2];
}
-void ztm(std::complex<double> **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6, C9 * c9) {
- complex<double> gie, gle, a1, a2, a3, a4, sum1, sum2, sum3, sum4;
- const complex<double> cc0(0.0, 0.0);
+void ztm(dcomplex **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6, C9 * c9) {
+ dcomplex gie, gle, a1, a2, a3, a4, sum1, sum2, sum3, sum4;
+ const dcomplex cc0 = 0.0 + 0.0 * I;
int ndi = c4->nsph * c4->nlim;
int i2 = 0;
for (int n2 = 1; n2 <= c4->nsph; n2++) { // GPU portable?
diff --git a/src/libnptm/lapack_calls.cpp b/src/libnptm/lapack_calls.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0b1dac5ac92ec6e876a560c0ab239e49a77314c2
--- /dev/null
+++ b/src/libnptm/lapack_calls.cpp
@@ -0,0 +1,45 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
+/*! \file lapack_calls.cpp
+ *
+ * \brief Implementation of the interface with LAPACK libraries.
+ */
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
+
+#ifdef USE_LAPACK
+#ifdef USE_MKL
+#include <mkl_lapacke.h>
+#else
+#include <lapacke.h>
+#endif
+
+#ifndef INCLUDE_LAPACK_CALLS_H_
+#include "../include/lapack_calls.h"
+#endif
+#endif
+
+#ifdef USE_LAPACK
+void zinvert(dcomplex **mat, lapack_int n, int &jer) {
+ jer = 0;
+ dcomplex *arr = &(mat[0][0]);
+ const dcomplex uim = 0.0 + 1.0 * I;
+
+#ifdef USE_MKL
+ MKL_Complex16 *arr2 = (MKL_Complex16 *) arr;
+#endif
+
+ lapack_int* IPIV = new lapack_int[n]();
+
+#ifdef USE_MKL
+ LAPACKE_zgetrf(LAPACK_ROW_MAJOR, n, n, arr2, n, IPIV);
+ LAPACKE_zgetri(LAPACK_ROW_MAJOR, n, arr2, n, IPIV);
+#else
+ LAPACKE_zgetrf(LAPACK_ROW_MAJOR, n, n, arr, n, IPIV);
+ LAPACKE_zgetri(LAPACK_ROW_MAJOR, n, arr, n, IPIV);
+#endif
+
+ delete[] IPIV;
+}
+#endif
diff --git a/src/libnptm/sph_subs.cpp b/src/libnptm/sph_subs.cpp
index 30e31bf9ec78e83a8e545a408d93c22fb6a267af..72f8cb1024912d78193dedbb3d71467e8e35de53 100644
--- a/src/libnptm/sph_subs.cpp
+++ b/src/libnptm/sph_subs.cpp
@@ -1,8 +1,12 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file sph_subs.cpp
*
* \brief C++ implementation of SPHERE subroutines.
*/
-#include <complex>
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
#ifndef INCLUDE_COMMONS_H_
#include "../include/Commons.h"
@@ -12,11 +16,9 @@
#include "../include/sph_subs.h"
#endif
-using namespace std;
-
void aps(double ****zpv, int li, int nsph, C1 *c1, double sqk, double *gaps) {
- complex<double> cc0 = complex<double>(0.0, 0.0);
- complex<double> summ, sume, suem, suee, sum;
+ dcomplex cc0 = 0.0 + 0.0 * I;
+ dcomplex summ, sume, suem, suee, sum;
double half_pi = acos(0.0);
double cofs = half_pi * 2.0 / sqk;
for (int i40 = 0; i40 < nsph; i40++) {
@@ -60,11 +62,11 @@ void aps(double ****zpv, int li, int nsph, C1 *c1, double sqk, double *gaps) {
}
}
}
- gaps[i40] = sum.real() * cofs;
+ gaps[i40] = real(sum) * cofs;
}
}
-void cbf(int n, std::complex<double> z, int &nm, std::complex<double> csj[]) {
+void cbf(int n, dcomplex z, int &nm, dcomplex *csj) {
/*
* FROM CSPHJY OF LIBRARY specfun
*
@@ -79,40 +81,40 @@ void cbf(int n, std::complex<double> z, int &nm, std::complex<double> csj[]) {
* point for backward recurrence
* ==========================================================
*/
- double zz = z.real() * z.real() + z.imag() * z.imag();
+ double zz = real(z) * real(z) + imag(z) * imag(z);
double a0 = sqrt(zz);
nm = n;
if (a0 < 1.0e-60) {
for (int k = 2; k <= n + 1; k++) {
csj[k - 1] = 0.0;
}
- csj[0] = complex<double>(1.0, 0.0);
+ csj[0] = 1.0 + 0.0 * I;
return;
}
- csj[0] = std::sin(z) / z;
+ csj[0] = csin(z) / z;
if (n == 0) {
return;
}
- csj[1] = (csj[0] -std::cos(z)) / z;
+ csj[1] = (csj[0] -ccos(z)) / z;
if (n == 1) {
return;
}
- complex<double> csa = csj[0];
- complex<double> csb = csj[1];
+ dcomplex csa = csj[0];
+ dcomplex csb = csj[1];
int m = msta1(a0, 200);
if (m < n) nm = m;
else m = msta2(a0, n, 15);
- complex<double> cf0 = 0.0;
- complex<double> cf1 = 1.0e-100;
- complex<double> cf, cs;
+ dcomplex cf0 = 0.0;
+ dcomplex cf1 = 1.0e-100;
+ dcomplex cf, cs;
for (int k = m; k >= 0; k--) {
cf = (2.0 * k + 3.0) * cf1 / z - cf0;
if (k <= nm) csj[k] = cf;
cf0 = cf1;
cf1 = cf;
}
- double abs_csa = abs(csa);
- double abs_csb = abs(csb);
+ double abs_csa = cabs(csa);
+ double abs_csb = cabs(csb);
if (abs_csa > abs_csb) cs = csa / cf;
else cs = csb / cf0;
for (int k = 0; k <= nm; k++) {
@@ -177,20 +179,20 @@ double cg1(int lmpml, int mu, int l, int m) {
return result;
}
-complex<double> dconjg(std::complex<double> z) {
- double zreal = z.real();
- double zimag = z.imag();
- return complex<double>(zreal, -zimag);
+dcomplex dconjg(dcomplex z) {
+ double zreal = real(z);
+ double zimag = imag(z);
+ return (zreal - zimag * I);
}
void diel(int npntmo, int ns, int i, int ic, double vk, C1 *c1, C2 *c2) {
const double dif = c1->rc[i - 1][ns] - c1->rc[i - 1][ns - 1];
const double half_step = 0.5 * dif / npntmo;
double rr = c1->rc[i - 1][ns - 1];
- const complex<double> delta = c2->dc0[ic] - c2->dc0[ic - 1];
+ const dcomplex delta = c2->dc0[ic] - c2->dc0[ic - 1];
const int kpnt = npntmo + npntmo;
c2->ris[kpnt] = c2->dc0[ic];
- c2->dlri[kpnt] = complex<double>(0.0, 0.0);
+ c2->dlri[kpnt] = 0.0 + 0.0 * I;
const int i90 = i - 1;
const int ns90 = ns - 1;
const int ic90 = ic - 1;
@@ -204,29 +206,29 @@ void diel(int npntmo, int ns, int i, int ic, double vk, C1 *c1, C2 *c2) {
void dme(
int li, int i, int npnt, int npntts, double vk, double exdc, double exri,
- C1 *c1, C2 *c2, int &jer, int &lcalc, std::complex<double> &arg
+ C1 *c1, C2 *c2, int &jer, int &lcalc, dcomplex &arg
) {
const int lipo = li + 1;
const int lipt = li + 2;
double *rfj = new double[lipt];
double *rfn = new double[lipt];
- complex<double> cfj[lipt], fbi[lipt], fb[lipt], fn[lipt];
- complex<double> rmf[li], drmf[li], ref[li], dref[li];
- complex<double> dfbi, dfb, dfn, ccna, ccnb, ccnc, ccnd;
- complex<double> y1, dy1, y2, dy2, arin, cri, uim;
+ dcomplex cfj[lipt], fbi[lipt], fb[lipt], fn[lipt];
+ dcomplex rmf[li], drmf[li], ref[li], dref[li];
+ dcomplex dfbi, dfb, dfn, ccna, ccnb, ccnc, ccnd;
+ dcomplex y1, dy1, y2, dy2, arin, cri;
+ const dcomplex uim = 1.0 * I;
jer = 0;
- uim = complex<double>(0.0, 1.0);
int nstp = npnt - 1;
int nstpts = npntts - 1;
double sz = vk * c1->ros[i - 1];
c2->vsz[i - 1] = sz;
double vkr1 = vk * c1->rc[i - 1][0];
int nsh = c1->nshl[i - 1];
- c2->vkt[i - 1] = std::sqrt(c2->dc0[0]);
+ c2->vkt[i - 1] = csqrt(c2->dc0[0]);
arg = vkr1 * c2->vkt[i - 1];
arin = arg;
bool goto32 = false;
- if (arg.imag() != 0.0) {
+ if (imag(arg) != 0.0) {
cbf(lipo, arg, lcalc, cfj);
if (lcalc < lipo) {
jer = 5;
@@ -238,7 +240,7 @@ void dme(
goto32 = true;
}
if (!goto32) {
- rbf(lipo, arg.real(), lcalc, rfj);
+ rbf(lipo, real(arg), lcalc, rfj);
if (lcalc < lipo) {
jer = 5;
delete[] rfj;
@@ -353,23 +355,23 @@ double envj(int n, double x) {
return result;
}
-void mmulc(std::complex<double> *vint, double **cmullr, double **cmul) {
- double sm2 = vint[0].real();
- double s24 = vint[1].real();
- double d24 = vint[1].imag();
- double sm4 = vint[5].real();
- double s23 = vint[8].real();
- double d32 = vint[8].imag();
- double s34 = vint[9].real();
- double d34 = vint[9].imag();
- double sm3 = vint[10].real();
- double s31 = vint[11].real();
- double d31 = vint[11].imag();
- double s21 = vint[12].real();
- double d12 = vint[12].imag();
- double s41 = vint[13].real();
- double d14 = vint[13].imag();
- double sm1 = vint[15].real();
+void mmulc(dcomplex *vint, double **cmullr, double **cmul) {
+ double sm2 = real(vint[0]);
+ double s24 = real(vint[1]);
+ double d24 = imag(vint[1]);
+ double sm4 = real(vint[5]);
+ double s23 = real(vint[8]);
+ double d32 = imag(vint[8]);
+ double s34 = real(vint[9]);
+ double d34 = imag(vint[9]);
+ double sm3 = real(vint[10]);
+ double s31 = real(vint[11]);
+ double d31 = imag(vint[11]);
+ double s21 = real(vint[12]);
+ double d12 = imag(vint[12]);
+ double s41 = real(vint[13]);
+ double d14 = imag(vint[13]);
+ double sm1 = real(vint[15]);
cmullr[0][0] = sm2;
cmullr[0][1] = sm3;
cmullr[0][2] = -s23;
@@ -462,7 +464,7 @@ int msta2(double x, int n, int mp) {
return result;
}
-void orunve( double *u1, double *u2, double *u3, int iorth, double torth) {
+void orunve(double *u1, double *u2, double *u3, int iorth, double torth) {
if (iorth <= 0) {
double cp = u1[0] * u2[0] + u1[1] * u2[1] + u1[2] * u2[2];
double abs_cp = cp;
@@ -481,7 +483,7 @@ void orunve( double *u1, double *u2, double *u3, int iorth, double torth) {
}
void pwma(
- double *up, double *un, std::complex<double> *ylm, int inpol, int lw,
+ double *up, double *un, dcomplex *ylm, int inpol, int lw,
int isq, C1 *c1
) {
const double four_pi = 8.0 * acos(0.0);
@@ -492,19 +494,19 @@ void pwma(
int nlwm = lw * (lw + 2);
int nlwmt = nlwm + nlwm;
const double sqrtwi = 1.0 / sqrt(2.0);
- const complex<double> uim(0.0, 1.0);
- complex<double> cm1 = 0.5 * complex<double>(up[0], up[1]);
- complex<double> cp1 = 0.5 * complex<double>(up[0], -up[1]);
+ const dcomplex uim = 1.0 * I;
+ dcomplex cm1 = 0.5 * (up[0] + up[1] * I);
+ dcomplex cp1 = 0.5 * (up[0] - up[1] * I);
double cz1 = up[2];
- complex<double> cm2 = 0.5 * complex<double>(un[0], un[1]);
- complex<double> cp2 = 0.5 * complex<double>(un[0], -un[1]);
+ dcomplex cm2 = 0.5 * (un[0] + un[1] * I);
+ dcomplex cp2 = 0.5 * (un[0] - un[1] * I);
double cz2 = un[2];
for (int l20 = 0; l20 < lw; l20++) {
int l = l20 + 1;
int lf = l + 1;
int lftl = lf * l;
double x = 1.0 * lftl;
- complex<double> cl = complex<double>(four_pi / sqrt(x), 0.0);
+ dcomplex cl = four_pi / sqrt(x) + 0.0 * I;
for (int powi = 1; powi <= l; powi++) cl *= uim;
int mv = l + lf;
int m = -lf;
@@ -536,8 +538,8 @@ void pwma(
if (inpol != 0) {
for (int k40 = 0; k40 < nlwm; k40++) {
int i = k40 + nlwm;
- complex<double> cc1 = sqrtwi * (c1->w[k40][ispo - 1] + uim * c1->w[k40][ispt - 1]);
- complex<double> cc2 = sqrtwi * (c1->w[k40][ispo - 1] - uim * c1->w[k40][ispt - 1]);
+ dcomplex cc1 = sqrtwi * (c1->w[k40][ispo - 1] + uim * c1->w[k40][ispt - 1]);
+ dcomplex cc2 = sqrtwi * (c1->w[k40][ispo - 1] - uim * c1->w[k40][ispt - 1]);
c1->w[k40][ispo - 1] = cc2;
c1->w[i][ispo - 1] = -cc2;
c1->w[k40][ispt - 1] = cc1;
@@ -560,11 +562,11 @@ void pwma(
}
void rabas(
- int inpol, int li, int nsph, C1 *c1, double **tqse, std::complex<double> **tqspe,
- double **tqss, std::complex<double> **tqsps
+ int inpol, int li, int nsph, C1 *c1, double **tqse, dcomplex **tqspe,
+ double **tqss, dcomplex **tqsps
) {
- complex<double> cc0 = complex<double>(0.0, 0.0);
- complex<double> uim = complex<double>(0.0, 1.0);
+ dcomplex cc0 = 0.0 + 0.0 * I;
+ dcomplex uim = 0.0 + 1.0 * I;
double two_pi = 4.0 * acos(0.0);
for (int i80 = 0; i80 < nsph; i80++) {
int i = i80 + 1;
@@ -580,19 +582,19 @@ void rabas(
for (int l70 = 0; l70 < li; l70++) {
int l = l70 + 1;
double fl = 1.0 * (l + l + 1);
- complex<double> rm = 1.0 / c1->rmi[l70][i80];
- double rmm = (rm * dconjg(rm)).real();
- complex<double> re = 1.0 / c1->rei[l70][i80];
- double rem = (re * dconjg(re)).real();
+ dcomplex rm = 1.0 / c1->rmi[l70][i80];
+ double rmm = real(rm * dconjg(rm));
+ dcomplex re = 1.0 / c1->rei[l70][i80];
+ double rem = real(re * dconjg(re));
if (inpol == 0) {
- complex<double> pce = ((rm + re) * uim) * fl;
- complex<double> pcs = ((rmm + rem) * fl) * uim;
+ dcomplex pce = ((rm + re) * uim) * fl;
+ dcomplex pcs = ((rmm + rem) * fl) * uim;
tqspe[0][i80] -= pce;
tqspe[1][i80] += pce;
tqsps[0][i80] -= pcs;
tqsps[1][i80] += pcs;
} else {
- double ce = (rm + re).real() * fl;
+ double ce = real(rm + re) * fl;
double cs = (rmm + rem) * fl;
tqse[0][i80] -= ce;
tqse[1][i80] += ce;
@@ -672,12 +674,11 @@ void rbf(int n, double x, int &nm, double sj[]) {
}
void rkc(
- int npntmo, double step, std::complex<double> dcc, double &x, int lpo,
- std::complex<double> &y1, std::complex<double> &y2, std::complex<double> &dy1,
- std::complex<double> &dy2
+ int npntmo, double step, dcomplex dcc, double &x, int lpo,
+ dcomplex &y1, dcomplex &y2, dcomplex &dy1, dcomplex &dy2
) {
- complex<double> cy1, cdy1, c11, cy23, yc2, c12, c13;
- complex<double> cy4, yy, c14, c21, c22, c23, c24;
+ dcomplex cy1, cdy1, c11, cy23, yc2, c12, c13;
+ dcomplex cy4, yy, c14, c21, c22, c23, c24;
double half_step = 0.5 * step;
double cl = 1.0 * lpo * (lpo - 1);
for (int ipnt60 = 0; ipnt60 < npntmo; ipnt60++) {
@@ -710,12 +711,11 @@ void rkc(
}
void rkt(
- int npntmo, double step, double &x, int lpo, std::complex<double> &y1,
- std::complex<double> &y2, std::complex<double> &dy1, std::complex<double> &dy2,
- C2 *c2
+ int npntmo, double step, double &x, int lpo, dcomplex &y1,
+ dcomplex &y2, dcomplex &dy1, dcomplex &dy2, C2 *c2
) {
- complex<double> cy1, cdy1, c11, cy23, cdy23, yc2, c12, c13;
- complex<double> cy4, cdy4, yy, c14, c21, c22, c23, c24;
+ dcomplex cy1, cdy1, c11, cy23, cdy23, yc2, c12, c13;
+ dcomplex cy4, cdy4, yy, c14, c21, c22, c23, c24;
double half_step = 0.5 * step;
double cl = 1.0 * lpo * (lpo - 1);
for (int ipnt60 = 0; ipnt60 < npntmo; ipnt60++) {
@@ -758,7 +758,7 @@ void rkt(
}
}
-void rnf(int n, double x, int &nm, double sy[]) {
+void rnf(int n, double x, int &nm, double *sy) {
/*
* FROM SPHJY OF LIBRARY specfun
*
@@ -804,7 +804,7 @@ void rnf(int n, double x, int &nm, double sy[]) {
void sphar(
double cosrth, double sinrth, double cosrph, double sinrph,
- int ll, std::complex<double> *ylm
+ int ll, dcomplex *ylm
) {
const int rmp_size = ll;
const int plegn_size = (ll + 1) * ll / 2 + ll + 1;
@@ -867,10 +867,10 @@ void sphar(
lmp = l0p + m;
save = pi4irs * plegn[lmp - 1];
lmy = l0y + m;
- ylm[lmy - 1] = save * complex<double>(cosrmp[m - 1], sinrmp[m - 1]);
+ ylm[lmy - 1] = save * (cosrmp[m - 1] + sinrmp[m - 1] * I);
if (m % 2 != 0) ylm[lmy - 1] *= -1.0;
lmy = l0y - m;
- ylm[lmy - 1] = save * complex<double>(cosrmp[m - 1], -sinrmp[m - 1]);
+ ylm[lmy - 1] = save * (cosrmp[m - 1] - sinrmp[m - 1] * I);
label45:
if (m >= l) goto label47;
m += 1;
@@ -881,33 +881,33 @@ void sphar(
goto label40;
}
-void sscr0(std::complex<double> &tfsas, int nsph, int lm, double vk, double exri, C1 *c1) {
- complex<double> sum21, rm, re, csam;
- const complex<double> cc0 = complex<double>(0.0, 0.0);
+void sscr0(dcomplex &tfsas, int nsph, int lm, double vk, double exri, C1 *c1) {
+ dcomplex sum21, rm, re, csam;
+ const dcomplex cc0 = 0.0 + 0.0 * I;
const double exdc = exri * exri;
double ccs = 4.0 * acos(0.0) / (vk * vk);
double cccs = ccs / exdc;
- csam = -(ccs / (exri * vk)) * complex<double>(0.0, 0.5);
+ csam = -(ccs / (exri * vk)) * (0.0 + 0.5 * I);
tfsas = cc0;
for (int i12 = 0; i12 < nsph; i12++) {
int i = i12 + 1;
int iogi = c1->iog[i12];
if (iogi >= i) {
double sums = 0.0;
- complex<double> sum21 = cc0;
+ dcomplex sum21 = cc0;
for (int l10 = 0; l10 < lm; l10++) {
int l = l10 + 1;
double fl = 1.0 + l + l;
rm = 1.0 / c1->rmi[l10][i12];
re = 1.0 / c1->rei[l10][i12];
- complex<double> rm_cnjg = dconjg(rm);
- complex<double> re_cnjg = dconjg(re);
- sums += (rm_cnjg * rm + re_cnjg * re).real() * fl;
+ dcomplex rm_cnjg = dconjg(rm);
+ dcomplex re_cnjg = dconjg(re);
+ sums += real(rm_cnjg * rm + re_cnjg * re) * fl;
sum21 += (rm + re) * fl;
}
sum21 *= -1.0;
double scasec = cccs * sums;
- double extsec = -cccs * sum21.real();
+ double extsec = -cccs * real(sum21);
double abssec = extsec - scasec;
c1->sscs[i12] = scasec;
c1->sexs[i12] = extsec;
@@ -923,10 +923,10 @@ void sscr0(std::complex<double> &tfsas, int nsph, int lm, double vk, double exri
}
void sscr2(int nsph, int lm, double vk, double exri, C1 *c1) {
- std::complex<double> s11, s21, s12, s22, rm, re, csam, cc;
- const complex<double> cc0(0.0, 0.0);
+ dcomplex s11, s21, s12, s22, rm, re, csam, cc;
+ const dcomplex cc0 = 0.0 + 0.0 * I;
double ccs = 1.0 / (vk * vk);
- csam = -(ccs / (exri * vk)) * complex<double>(0.0, 0.5);
+ csam = -(ccs / (exri * vk)) * (0.0 + 0.5 * I);
const double pigfsq = 64.0 * acos(0.0) * acos(0.0);
double cfsq = 4.0 / (pigfsq * ccs * ccs);
int nlmm = lm * (lm + 2);
@@ -966,7 +966,7 @@ void sscr2(int nsph, int lm, double vk, double exri, C1 *c1) {
int j = 0;
for (int ipo1 = 0; ipo1 < 2; ipo1++) {
for (int jpo1 = 0; jpo1 < 2; jpo1++) {
- complex<double> cc = dconjg(c1->sas[i24][jpo1][ipo1]);
+ dcomplex cc = dconjg(c1->sas[i24][jpo1][ipo1]);
for (int ipo2 = 0; ipo2 < 2; ipo2++) {
for (int jpo2 = 0; jpo2 < 2; jpo2++) {
c1->vints[i24][j++] = c1->sas[i24][jpo2][ipo2] * cc * cfsq;
@@ -1111,9 +1111,9 @@ void wmamp(
double *un, C1 *c1
) {
const int ylm_size = (lm + 1) * (lm + 1) + 1;
- complex<double> *ylm = new complex<double>[ylm_size];
+ dcomplex *ylm = new dcomplex[ylm_size];
const int nlmp = lm * (lm + 2) + 2;
- ylm[nlmp - 1] = complex<double>(0.0, 0.0);
+ ylm[nlmp - 1] = 0.0 + 0.0 * I;
if (idot != 0) {
if (idot != 1) {
for (int n40 = 0; n40 < nsph; n40++) {
@@ -1140,9 +1140,9 @@ void wmasp(
int nsph, double *argi, double *args, C1 *c1
) {
const int ylm_size = (lm + 1) * (lm + 1) + 1;
- complex<double> *ylm = new complex<double>[ylm_size];
+ dcomplex *ylm = new dcomplex[ylm_size];
const int nlmp = lm * (lm + 2) + 2;
- ylm[nlmp - 1] = complex<double>(0.0, 0.0);
+ ylm[nlmp - 1] = 0.0 + 0.0 * I;
if (idot != 0) {
if (idot != 1) {
for (int n40 = 0; n40 < nsph; n40++) {
diff --git a/src/libnptm/tfrfme.cpp b/src/libnptm/tfrfme.cpp
index b3478d3cc34850289b6a583a6db78fbd0ae0f663..a1a8e1aecf19841471964057b2d3442711638e82 100644
--- a/src/libnptm/tfrfme.cpp
+++ b/src/libnptm/tfrfme.cpp
@@ -1,14 +1,19 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file tfrfme.cpp
*
* \brief Implementation of the trapping calculation objects.
*/
-#include <complex>
#include <exception>
#include <fstream>
#include <hdf5.h>
#include <string>
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
+
#ifndef INCLUDE_ERRORS_H_
#include "../include/errors.h"
#endif
@@ -32,7 +37,7 @@ Swap1::Swap1(int lm, int _nkv) {
nkv = _nkv;
nlmmt = 2 * lm * (lm + 2);
const int size = nkv * nkv * nlmmt;
- wk = new complex<double>[size]();
+ wk = new dcomplex[size]();
last_index = 0;
}
@@ -57,8 +62,8 @@ Swap1* Swap1::from_hdf5(string file_name) {
double *elements;
string str_type;
int _nlmmt, _nkv, lm, num_elements, index;
- complex<double> value;
- complex<double> *_wk = NULL;
+ dcomplex value;
+ dcomplex *_wk = NULL;
if (status == 0) {
status = hdf_file->read("NLMMT", "INT32", &_nlmmt);
status = hdf_file->read("NKV", "INT32", &_nkv);
@@ -71,7 +76,7 @@ Swap1* Swap1::from_hdf5(string file_name) {
status = hdf_file->read("WK", str_type, elements);
for (int wi = 0; wi < num_elements / 2; wi++) {
index = 2 * wi;
- value = complex<double>(elements[index], elements[index + 1]);
+ value = elements[index] + elements[index + 1] * I;
_wk[wi] = value;
} // wi loop
delete[] elements;
@@ -84,7 +89,7 @@ Swap1* Swap1::from_hdf5(string file_name) {
Swap1* Swap1::from_legacy(string file_name) {
fstream input;
Swap1 *instance = NULL;
- complex<double> *_wk = NULL;
+ dcomplex *_wk = NULL;
int _nlmmt, _nkv, lm;
double rval, ival;
input.open(file_name.c_str(), ios::in | ios::binary);
@@ -98,7 +103,7 @@ Swap1* Swap1::from_legacy(string file_name) {
for (int j = 0; j < num_elements; j++) {
input.read(reinterpret_cast<char *>(&rval), sizeof(double));
input.read(reinterpret_cast<char *>(&ival), sizeof(double));
- _wk[j] = complex<double>(rval, ival);
+ _wk[j] = rval + ival * I;
}
input.close();
} else {
@@ -109,7 +114,7 @@ Swap1* Swap1::from_legacy(string file_name) {
long Swap1::get_memory_requirement(int lm, int _nkv) {
long size = (long)(3 * sizeof(int));
- size += (long)(sizeof(complex<double>) * 2 * lm * (lm + 2) * _nkv * _nkv);
+ size += (long)(sizeof(dcomplex) * 2 * lm * (lm + 2) * _nkv * _nkv);
return size;
}
@@ -142,8 +147,8 @@ void Swap1::write_hdf5(string file_name) {
rec_type_list.append(str_type);
double *ptr_elements = new double[num_elements]();
for (int wi = 0; wi < num_elements / 2; wi++) {
- ptr_elements[2 * wi] = wk[wi].real();
- ptr_elements[2 * wi + 1] = wk[wi].imag();
+ ptr_elements[2 * wi] = real(wk[wi]);
+ ptr_elements[2 * wi + 1] = imag(wk[wi]);
}
rec_ptr_list.append(ptr_elements);
@@ -173,8 +178,8 @@ void Swap1::write_legacy(string file_name) {
output.write(reinterpret_cast<char *>(&nlmmt), sizeof(int));
output.write(reinterpret_cast<char *>(&nkv), sizeof(int));
for (int j = 0; j < num_elements; j++) {
- rval = wk[j].real();
- ival = wk[j].imag();
+ rval = real(wk[j]);
+ ival = imag(wk[j]);
output.write(reinterpret_cast<char *>(&rval), sizeof(double));
output.write(reinterpret_cast<char *>(&ival), sizeof(double));
}
@@ -591,8 +596,8 @@ TFRFME::TFRFME(int _lmode, int _lm, int _nkv, int _nxv, int _nyv, int _nzv) {
zv = new double[nzv]();
nlmmt = lm * (lm + 2) * 2;
nrvc = nxv * nyv * nzv;
- wsum = new complex<double>*[nlmmt];
- for (int wi = 0; wi < nlmmt; wi++) wsum[wi] = new complex<double>[nrvc]();
+ wsum = new dcomplex*[nlmmt];
+ for (int wi = 0; wi < nlmmt; wi++) wsum[wi] = new dcomplex[nrvc]();
}
TFRFME::~TFRFME() {
@@ -624,8 +629,8 @@ TFRFME* TFRFME::from_hdf5(string file_name) {
double *elements;
string str_type;
int _nlmmt, _nrvc, num_elements;
- complex<double> value;
- complex<double> **_wsum = NULL;
+ dcomplex value;
+ dcomplex **_wsum = NULL;
if (status == 0) {
int lmode, lm, nkv, nxv, nyv, nzv;
double vk, exri, an, ff, tra, spd, frsh, exril;
@@ -668,7 +673,7 @@ TFRFME* TFRFME::from_hdf5(string file_name) {
int index = 0;
for (int wj = 0; wj < _nrvc; wj++) {
for (int wi = 0; wi < _nlmmt; wi++) {
- value = complex<double>(elements[index], elements[index + 1]);
+ value = elements[index] + elements[index + 1] * I;
_wsum[wi][wj] = value;
index += 2;
} // wi loop
@@ -683,7 +688,7 @@ TFRFME* TFRFME::from_hdf5(string file_name) {
TFRFME* TFRFME::from_legacy(string file_name) {
fstream input;
TFRFME *instance = NULL;
- complex<double> **_wsum = NULL;
+ dcomplex **_wsum = NULL;
double *coord_vec = NULL;
input.open(file_name.c_str(), ios::in | ios::binary);
if (input.is_open()) {
@@ -735,7 +740,7 @@ TFRFME* TFRFME::from_legacy(string file_name) {
for (int wi = 0; wi < _nlmmt; wi++) {
input.read(reinterpret_cast<char *>(&rval), sizeof(double));
input.read(reinterpret_cast<char *>(&ival), sizeof(double));
- complex<double> value(rval, ival);
+ dcomplex value = rval + ival * I;
_wsum[wi][wj] = value;
} // wi loop
} // wj loop
@@ -754,7 +759,7 @@ long TFRFME::get_memory_requirement(
int _nrvc = _nxv * _nyv * _nzv;
long size = (long)(8 * sizeof(int) + 8 * sizeof(double));
size += (long)((_nxv + _nyv + _nzv) * sizeof(double));
- size += (long)(_nlmmt * _nrvc * sizeof(complex<double>));
+ size += (long)(_nlmmt * _nrvc * sizeof(dcomplex));
return size;
}
@@ -877,8 +882,8 @@ void TFRFME::write_hdf5(string file_name) {
int index = 0;
for (int wj = 0; wj < nrvc; wj++) {
for (int wi = 0; wi < nlmmt; wi++) {
- ptr_elements[index++] = wsum[wi][wj].real();
- ptr_elements[index++] = wsum[wi][wj].imag();
+ ptr_elements[index++] = real(wsum[wi][wj]);
+ ptr_elements[index++] = imag(wsum[wi][wj]);
} // wi loop
} // wj loop
rec_ptr_list.append(ptr_elements);
@@ -926,8 +931,8 @@ void TFRFME::write_legacy(string file_name) {
output.write(reinterpret_cast<char *>(&(zv[zi])), sizeof(double));
for (int wj = 0; wj < nrvc; wj++) {
for (int wi = 0; wi < nlmmt; wi++) {
- double rval = wsum[wi][wj].real();
- double ival = wsum[wi][wj].imag();
+ double rval = real(wsum[wi][wj]);
+ double ival = imag(wsum[wi][wj]);
output.write(reinterpret_cast<char *>(&rval), sizeof(double));
output.write(reinterpret_cast<char *>(&ival), sizeof(double));
} // wi loop
diff --git a/src/libnptm/tra_subs.cpp b/src/libnptm/tra_subs.cpp
index 47c2b28edfc1fb784f06712f8f70baff9830b5ad..2591b8c335a245a64cd41ae26ad50e6720810ace 100644
--- a/src/libnptm/tra_subs.cpp
+++ b/src/libnptm/tra_subs.cpp
@@ -1,11 +1,16 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file tra_subs.cpp
*
* \brief C++ implementation of TRAPPING subroutines.
*/
#include <cmath>
-#include <complex>
#include <fstream>
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
+
#ifndef INCLUDE_COMMONS_H_
#include "../include/Commons.h"
#endif
@@ -22,12 +27,7 @@
#include "../include/tra_subs.h"
#endif
-using namespace std;
-
-void camp(
- std::complex<double> *ac, std::complex<double> **am0m, std::complex<double> *ws,
- CIL *cil
-) {
+void camp(dcomplex *ac, dcomplex **am0m, dcomplex *ws, CIL *cil) {
for (int j = 0; j < cil->nlemt; j++) {
for (int i = 0; i < cil->nlemt; i++) {
ac[j] += (am0m[j][i] * ws[i]);
@@ -35,13 +35,10 @@ void camp(
} // j loop
}
-void czamp(
- std::complex<double> *ac, std::complex<double> **amd, int **indam,
- std::complex<double> *ws, CIL *cil
-) {
- const complex<double> cc0(0.0, 0.0);
- const complex<double> uim(0.0, 1.0);
- complex<double> summ, sume;
+void czamp(dcomplex *ac, dcomplex **amd, int **indam, dcomplex *ws, CIL *cil) {
+ const dcomplex cc0 = 0.0 + 0.0 * I;
+ const dcomplex uim = 0.0 + 1.0 * I;
+ dcomplex summ, sume;
for (int im20 = 1; im20 <= cil->mxim; im20++) {
int m = im20 - cil->mxmpo;
int abs_m = (m < 0) ? -m : m;
@@ -65,13 +62,13 @@ void czamp(
}
void ffrf(
- double ****zpv, std::complex<double> *ac, std::complex<double> *ws, double *fffe,
+ double ****zpv, dcomplex *ac, dcomplex *ws, double *fffe,
double *fffs, CIL *cil, CCR *ccr
) {
- const complex<double> cc0(0.0, 0.0);
- const complex<double> uim(0.0, 1.0);
- complex<double> uimmp, summ, sume, suem, suee;
- complex<double> *gap = new complex<double>[3]();
+ const dcomplex cc0 = 0.0 + 0.0 * I;
+ const dcomplex uim = 0.0 + 1.0 * I;
+ dcomplex uimmp, summ, sume, suem, suee;
+ dcomplex *gap = new dcomplex[3]();
for (int imu50 = 1; imu50 <= 3; imu50++) {
int mu = imu50 - 2;
@@ -123,9 +120,9 @@ void ffrf(
sume = -gap[0] * ccr->cimu;
suee = -gap[1] * ccr->cof;
suem = -gap[2] * ccr->cimu;
- fffe[0] = (sume - suem).real();
- fffe[1] = ((sume + suem) * uim).real();
- fffe[2] = suee.real();
+ fffe[0] = real(sume - suem);
+ fffe[1] = real((sume + suem) * uim);
+ fffe[2] = real(suee);
for (int imu90 = 1; imu90 <= 3; imu90++) {
int mu = imu90 - 2;
@@ -177,25 +174,22 @@ void ffrf(
sume = gap[0] * ccr->cimu;
suee = gap[1] * ccr->cof;
suem = gap[2] * ccr->cimu;
- fffs[0] = (sume - suem).real();
- fffs[1] = ((sume + suem) * uim).real();
- fffs[2] = suee.real();
+ fffs[0] = real(sume - suem);
+ fffs[1] = real((sume + suem) * uim);
+ fffs[2] = real(suee);
delete[] gap;
}
-void ffrt(
- std::complex<double> *ac, std::complex<double> *ws, double *ffte, double *ffts,
- CIL *cil
-) {
- const complex<double> cc0(0.0, 0.0);
- const complex<double> uim(0.0, 1.0);
+void ffrt(dcomplex *ac, dcomplex *ws, double *ffte, double *ffts, CIL *cil) {
+ const dcomplex cc0 = 0.0 + 0.0 * I;
+ const dcomplex uim = 0.0 + 1.0 * I;
const double sq2i = 1.0 / sqrt(2.0);
- const complex<double> sq2iti = uim * sq2i;
- complex<double> aca, acw;
- complex<double> *ctqce, *ctqcs;
+ const dcomplex sq2iti = uim * sq2i;
+ dcomplex aca, acw;
+ dcomplex *ctqce, *ctqcs;
- ctqce = new complex<double>[3]();
- ctqcs = new complex<double>[3]();
+ ctqce = new dcomplex[3]();
+ ctqcs = new dcomplex[3]();
for (int l60 = 1; l60 <= cil->le; l60++) {
int lpo = l60 + 1;
int il = l60 * lpo;
@@ -235,28 +229,28 @@ void ffrt(
}
} // im60 loop
} // l60 loop
- ffte[0] = (ctqce[0] - ctqce[2]).real() * sq2i;
- ffte[1] = (sq2iti * (ctqce[0] + ctqce[2])).real();
- ffte[2] = ctqce[1].real();
- ffts[0] = -sq2i * (ctqcs[0] - ctqcs[2]).real();
- ffts[1] = -1.0 * (sq2iti * (ctqcs[0] + ctqcs[2])).real();
- ffts[2] = -1.0 * ctqcs[1].real();
+ ffte[0] = real(ctqce[0] - ctqce[2]) * sq2i;
+ ffte[1] = real(sq2iti * (ctqce[0] + ctqce[2]));
+ ffte[2] = real(ctqce[1]);
+ ffts[0] = -sq2i * real(ctqcs[0] - ctqcs[2]);
+ ffts[1] = -1.0 * real(sq2iti * (ctqcs[0] + ctqcs[2]));
+ ffts[2] = -1.0 * real(ctqcs[1]);
delete[] ctqce;
delete[] ctqcs;
}
-complex<double> *frfmer(
+dcomplex *frfmer(
int nkv, double vkm, double vknmx, double apfafa, double tra,
double spd, double rir, double ftcn, int le, int lmode, double pmf,
Swap1 *tt1, Swap2 *tt2
) {
const int nlemt = le * (le + 2) * 2;
- const complex<double> cc0(0.0, 0.0);
+ const dcomplex cc0 = 0.0 + 0.0 * I;
double sq = vkm * vkm;
double *_vkv = tt2->get_vector();
double **_vkzm = tt2->get_matrix();
- complex<double> *wk = new complex<double>[nlemt]();
+ dcomplex *wk = new dcomplex[nlemt]();
for (int jy90 = 0; jy90 < nkv; jy90++) {
double vky = _vkv[jy90];
double sqy = vky * vky;
@@ -282,22 +276,22 @@ complex<double> *frfmer(
return wk;
}
-void pwmalp(std::complex<double> **w, double *up, double *un, std::complex<double> *ylm, int lw) {
- complex<double> cp1, cm1, cp2, cm2, cl;
- const complex<double> uim(0.0, 1.0);
+void pwmalp(dcomplex **w, double *up, double *un, dcomplex *ylm, int lw) {
+ dcomplex cp1, cm1, cp2, cm2, cl;
+ const dcomplex uim = 0.0 + 1.0 * I;
const double four_pi = acos(0.0) * 8.0;
const int nlwm = lw * (lw + 2);
- cm1 = 0.5 * complex<double>(up[0], up[1]);
- cp1 = 0.5 * complex<double>(up[0], -up[1]);
+ cm1 = 0.5 * (up[0] + up[1] * I);
+ cp1 = 0.5 * (up[0] - up[1] * I);
double cz1 = up[2];
- cm2 = 0.5 * complex<double>(un[0], un[1]);
- cp2 = 0.5 * complex<double>(un[0], -un[1]);
+ cm2 = 0.5 * (un[0] + un[1] * I);
+ cp2 = 0.5 * (un[0] - un[1] * I);
double cz2 =un[2];
for (int l20 = 1; l20 <= lw; l20++) {
int lf = l20 + 1;
int lftl = lf * l20;
double x = 1.0 * lftl;
- complex<double> cl = (four_pi / sqrt(x)) * std::pow(uim, 1.0 * l20);
+ dcomplex cl = (four_pi / sqrt(x)) * cpow(uim, 1.0 * l20);
int mv = l20 + lf;
int m = -lf;
for (int mf20 = 1; mf20 <= mv; mf20++) {
@@ -319,10 +313,7 @@ void pwmalp(std::complex<double> **w, double *up, double *un, std::complex<doubl
} // k30 loop
}
-void samp(
- std::complex<double> *ac, std::complex<double> *tmsm, std::complex<double> *tmse,
- std::complex<double> *ws, CIL *cil
-) {
+void samp(dcomplex *ac, dcomplex *tmsm, dcomplex *tmse, dcomplex *ws, CIL *cil) {
int i = 0;
for (int l20 = 0; l20 < cil->le; l20++) {
int l = l20 + 1;
@@ -336,13 +327,10 @@ void samp(
} // l20 loop
}
-void sampoa(
- std::complex<double> *ac, std::complex<double> **tms, std::complex<double> *ws,
- CIL *cil
-) {
- complex<double> **tm = new complex<double>*[2];
- tm[0] = new complex<double>[2]();
- tm[1] = new complex<double>[2]();
+void sampoa(dcomplex *ac, dcomplex **tms, dcomplex *ws, CIL *cil) {
+ dcomplex **tm = new dcomplex*[2];
+ tm[0] = new dcomplex[2]();
+ tm[1] = new dcomplex[2]();
int i = 0;
for (int l20 = 0; l20 < cil->le; l20++) {
tm[0][0] = tms[l20][0];
@@ -364,23 +352,23 @@ void sampoa(
}
void wamff(
- std::complex<double> *wk, double x, double y, double z, int lm, double apfafa,
+ dcomplex *wk, double x, double y, double z, int lm, double apfafa,
double tra, double spd, double rir, double ftcn, int lmode, double pmf
) {
const int nlmm = lm * (lm + 2);
const int nlmmt = 2 * nlmm;
const int nlmp = nlmm + 2;
- complex<double> **w, *ylm;
- const complex<double> cc0(0.0, 0.0);
- const complex<double> uim(0.0, 1.0);
- complex<double> cfam, cf1, cf2;
+ dcomplex **w, *ylm;
+ const dcomplex cc0 = 0.0 + 0.0 * I;
+ const dcomplex uim = 0.0 + 1.0 * I;
+ dcomplex cfam, cf1, cf2;
double rho, cph, sph, cth, sth, r;
double ftc1, ftc2;
double *up = new double[3];
double *un = new double[3];
- w = new complex<double>*[nlmmt];
- for (int wi = 0; wi < nlmmt; wi++) w[wi] = new complex<double>[2]();
- ylm = new complex<double>[nlmp]();
+ w = new dcomplex*[nlmmt];
+ for (int wi = 0; wi < nlmmt; wi++) w[wi] = new dcomplex[2]();
+ ylm = new dcomplex[nlmp]();
bool onx = (y == 0.0);
bool ony = (x == 0.0);
bool onz = (onx && ony);
@@ -446,7 +434,7 @@ void wamff(
double sthl = sth * rir;
double cthl = sqrt(1.0 - sthl * sthl);
double arg = spd * (z - (r / rir) * cthl);
- cfam = (tra * std::exp(-apfa * apfa) / sqrt(cthl)) * std::exp(uim * arg);
+ cfam = (tra * cexp(-apfa * apfa) / sqrt(cthl)) * cexp(uim * arg);
if (lmode == 0) {
if (sth == 0.0) { // label 45
ftc1 = ftcn;
@@ -478,7 +466,7 @@ void wamff(
skip62 = false;
}
} else { // label 57
- double fam = tra * std::exp(-apfa * apfa) / sqrt(cth);
+ double fam = tra * exp(-apfa * apfa) / sqrt(cth);
if (lmode == 0) {
double f1 = cph * fam;
double f2 = -sph * fam;
diff --git a/src/libnptm/types.cpp b/src/libnptm/types.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a9fd71be5682403b831ca13926fb6eff88ebac69
--- /dev/null
+++ b/src/libnptm/types.cpp
@@ -0,0 +1,13 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
+/*! \file types.cpp
+ *
+ * \brief Implementation of the functions connected with types.
+ */
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
+
+double real(dcomplex z) { return __real__ z; }
+
+double imag(dcomplex z) { return __imag__ z; }
diff --git a/src/make.inc b/src/make.inc
index dcdd896f0619de451fcf37eaeb62ec664f21c8e6..2ba3d6b135da8c1627da9b0e5504c6c266cc8777 100644
--- a/src/make.inc
+++ b/src/make.inc
@@ -38,10 +38,48 @@ ifndef HDF5_INCLUDE
override HDF5_INCLUDE=/usr/include/hdf5/serial
endif
+# define (outside) USE_LAPACK for lapacke support, LAPACK_ILP64 for ilp64 interface, MKL_ILP64 the same if using MKL implementation
+ifdef USE_LAPACK
+ifndef LAPACK_ILP64
+override LAPACK_ILP64=1
+endif
+# define (outside) USE_MKL to use the MKL implementation of lapacke
+ifdef USE_MKL
+ifndef MKL_ILP64
+override MKL_ILP64=1
+endif
+ifndef LAPACK_INCLUDE
+# this is for the MKL implementation
+override LAPACK_INCLUDE=$(MKLROOT)/include
+endif
+ifndef LAPACK_LDFLAGS
+# this is for the MKL implementation
+override LAPACK_LDFLAGS=-L$(MKLROOT)/lib -Wl,--no-as-needed -lmkl_intel_ilp64 -lmkl_gnu_thread -lmkl_core -lgomp -lpthread -lm -ldl
+endif
+# the next else refers to USE_MKL
+else
+ifndef LAPACK_INCLUDE
+# this is for standard "vanilla" lapacke64
+override LAPACK_INCLUDE=/usr/include
+endif
+ifndef LAPACK_LDFLAGS
+# this is for standard "vanilla" lapacke64
+override LAPACK_LDFLAGS=-llapacke64
+endif
+# the next endif is for USE_MKL
+endif
+#the next endif is for USE_LAPACK
+endif
+
# CXXFLAGS defines the default compilation options for the C++ compiler
ifndef CXXFLAGS
override CXXFLAGS=-O3 -ggdb -pg -coverage -I$(HDF5_INCLUDE)
-#override CXXFLAGS=-O3 -I$(HDF5_INCLUDE)
+ifdef USE_LAPACK
+override CXXFLAGS+= -DUSE_LAPACK -DLAPACK_ILP64
+ifdef USE_MKL
+override CXXFLAGS+= -DMKL_ILP64 -DUSE_MKL -I$(MKLROOT)/include
+endif
+endif
endif
# HDF5_LIB defines the default path to the HDF5 libraries to use
@@ -50,9 +88,11 @@ ifndef CXXLDFLAGS
ifndef HDF5_LIB
override HDF5_LIB=/usr/lib/x86_64-linux-gnu/hdf5/serial
endif
-override CXXLDFLAGS=-L/usr/lib64 -L$(HDF5_LIB) -lhdf5 $(LDFLAGS)
-#else
-#override CXXLDFLAGS=-L/usr/lib64 -L$(HDF5_LIB) -lhdf5 $(CXXLDFLAGS)
+override CXXLDFLAGS=-L/usr/lib64 -L$(HDF5_LIB) -lhdf5
+ifdef USE_LAPACK
+override CXXLDFLAGS+= $(LAPACK_LDFLAGS)
+endif
+override CXXLDFLAGS+= $(LDFLAGS)
endif
#SOFLAGS defines the additional flags for the c++ compiler to create a shared object file
diff --git a/src/scripts/README.md b/src/scripts/README.md
index b865eacaa0946125f0a30195fba197fc3f6c385d..59a757141c8c6627065b4b809e774a0bc9f079fd 100644
--- a/src/scripts/README.md
+++ b/src/scripts/README.md
@@ -1,3 +1,5 @@
+Distributed under the terms of *GNU GPL-3.0-or-later*
+
# Folder instructions
This directory contains scripts and tools to evaluate the code functionality.
diff --git a/src/sphere/np_sphere.cpp b/src/sphere/np_sphere.cpp
index e3033dae6fa36f373e17c1fda92f802537ce3b53..7f381e41ee28b645b53376828ddf33c4c68b7a88 100644
--- a/src/sphere/np_sphere.cpp
+++ b/src/sphere/np_sphere.cpp
@@ -1,3 +1,5 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file np_sphere.cpp
*
* \brief Single sphere scattering problem handler.
@@ -15,10 +17,13 @@
* while the advanced mode implements the possibility to change input and
* output options, without having to modify the code.
*/
-#include <complex>
#include <cstdio>
#include <string>
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
+
#ifndef INCLUDE_CONFIGURATION_H_
#include "../include/Configuration.h"
#endif
diff --git a/src/sphere/sphere.cpp b/src/sphere/sphere.cpp
index 7016d8361eb8d61fef8437f504516352f0444948..b4bf9652d3e0a9401e076f750b92615073750b97 100644
--- a/src/sphere/sphere.cpp
+++ b/src/sphere/sphere.cpp
@@ -1,13 +1,18 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file sphere.cpp
*
* \brief Implementation of the single sphere calculation.
*/
-#include <complex>
#include <cstdio>
#include <exception>
#include <fstream>
#include <string>
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
+
#ifndef INCLUDE_ERRORS_H_
#include "../include/errors.h"
#endif
@@ -37,7 +42,7 @@ using namespace std;
* \param output_path: `string` Directory to write the output files in.
*/
void sphere(string config_file, string data_file, string output_path) {
- complex<double> arg, s0, tfsas;
+ dcomplex arg, s0, tfsas;
double th, ph;
printf("INFO: making legacy configuration...\n");
ScattererConfiguration *sconf = NULL;
@@ -57,7 +62,7 @@ void sphere(string config_file, string data_file, string output_path) {
} catch(const OpenConfigurationFileException &ex) {
printf("\nERROR: failed to open geometry configuration file.\n");
printf("FILE: %s\n", ex.what());
- if (sconf) delete sconf;
+ if (sconf != NULL) delete sconf;
exit(1);
}
if (sconf->number_of_spheres == gconf->number_of_spheres) {
@@ -83,21 +88,21 @@ void sphere(string config_file, string data_file, string output_path) {
cmul[i] = new double[4];
cmullr[i] = new double[4];
}
- complex<double> **tqspe, **tqsps;
+ dcomplex **tqspe, **tqsps;
double **tqse, **tqss;
tqse = new double*[2];
tqss = new double*[2];
- tqspe = new std::complex<double>*[2];
- tqsps = new std::complex<double>*[2];
+ tqspe = new dcomplex*[2];
+ tqsps = new dcomplex*[2];
for (int ti = 0; ti < 2; ti++) {
tqse[ti] = new double[2]();
tqss[ti] = new double[2]();
- tqspe[ti] = new std::complex<double>[2]();
- tqsps[ti] = new std::complex<double>[2]();
+ tqspe[ti] = new dcomplex[2]();
+ tqsps[ti] = new dcomplex[2]();
}
double frx = 0.0, fry = 0.0, frz = 0.0;
double cfmp, cfsp, sfmp, sfsp;
- complex<double> *vint = new complex<double>[16];
+ dcomplex *vint = new dcomplex[16];
int jw;
int nsph = gconf->number_of_spheres;
C1 *c1 = new C1(nsph, gconf->l_max, sconf->nshl_vec, sconf->iog_vec);
@@ -244,8 +249,8 @@ void sphere(string config_file, string data_file, string output_path) {
fprintf(output, " \n");
}
for (int jxi488 = 0; jxi488 < sconf->number_of_scales; jxi488++) {
- printf("INFO: running scale iteration...\n");
int jxi = jxi488 + 1;
+ printf("INFO: running scale iteration %d of %d...\n", jxi, sconf->number_of_scales);
fprintf(output, "========== JXI =%3d ====================\n", jxi);
double xi = sconf->scale_vec[jxi488];
if (sconf->idfc >= 0) {
@@ -289,7 +294,7 @@ void sphere(string config_file, string data_file, string output_path) {
);
if (jer != 0) {
fprintf(output, " STOP IN DME\n");
- fprintf(output, " AT %1d LCALC=%3d TOO SMALL WITH ARG=%15.7lE+i(%15.7lE)\n", jer, lcalc, arg.real(), arg.imag());
+ fprintf(output, " AT %1d LCALC=%3d TOO SMALL WITH ARG=%15.7lE+i(%15.7lE)\n", jer, lcalc, real(arg), imag(arg));
tppoan.close();
fclose(output);
return;
@@ -323,8 +328,8 @@ void sphere(string config_file, string data_file, string output_path) {
output,
" SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE%15.7lE\n",
c2->vsz[i170],
- c2->vkt[i170].real(),
- c2->vkt[i170].imag()
+ real(c2->vkt[i170]),
+ imag(c2->vkt[i170])
);
}
fprintf(output, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n");
@@ -341,12 +346,12 @@ void sphere(string config_file, string data_file, string output_path) {
c1->sqscs[i170], c1->sqabs[i170],
c1->sqexs[i170]
);
- fprintf(output, " FSAS=%15.7lE%15.7lE\n", c1->fsas[i170].real(), c1->fsas[i170].imag());
+ fprintf(output, " FSAS=%15.7lE%15.7lE\n", real(c1->fsas[i170]), imag(c1->fsas[i170]));
double csch = 2.0 * vk * sqsfi / c1->gcsv[i170];
s0 = c1->fsas[i170] * exri;
- double qschu = csch * s0.imag();
- double pschu = csch * s0.real();
- double s0mag = cs0 * abs(s0);
+ double qschu = csch * imag(s0);
+ double pschu = csch * real(s0);
+ double s0mag = cs0 * cabs(s0);
fprintf(
output,
" QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",
@@ -359,33 +364,33 @@ void sphere(string config_file, string data_file, string output_path) {
fprintf(output, " IPO=%2d, TQEk=%15.7lE, TQSk=%15.7lE\n", 2, tqse[1][i170], tqss[1][i170]);
tppoan.write(reinterpret_cast<char *>(&(tqse[0][i170])), sizeof(double));
tppoan.write(reinterpret_cast<char *>(&(tqss[0][i170])), sizeof(double));
- double val = tqspe[0][i170].real();
+ double val = real(tqspe[0][i170]);
tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
- val = tqspe[0][i170].imag();
+ val = imag(tqspe[0][i170]);
tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
- val = tqsps[0][i170].real();
+ val = real(tqsps[0][i170]);
tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
- val = tqsps[0][i170].imag();
+ val = imag(tqsps[0][i170]);
tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
tppoan.write(reinterpret_cast<char *>(&(tqse[1][i170])), sizeof(double));
tppoan.write(reinterpret_cast<char *>(&(tqss[1][i170])), sizeof(double));
- val = tqspe[1][i170].real();
+ val = real(tqspe[1][i170]);
tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
- val = tqspe[1][i170].imag();
+ val = imag(tqspe[1][i170]);
tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
- val = tqsps[1][i170].real();
+ val = real(tqsps[1][i170]);
tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
- val = tqsps[1][i170].imag();
+ val = imag(tqsps[1][i170]);
tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
} // End if iog[i170] >= i
} // i170 loop
if (nsph != 1) {
- fprintf(output, " FSAT=(%15.7lE,%15.7lE)\n", tfsas.real(), tfsas.imag());
+ fprintf(output, " FSAT=(%15.7lE,%15.7lE)\n", real(tfsas), imag(tfsas));
double csch = 2.0 * vk * sqsfi / gcs;
s0 = tfsas * exri;
- double qschu = csch * s0.imag();
- double pschu = csch * s0.real();
- double s0mag = cs0 * abs(s0);
+ double qschu = csch * imag(s0);
+ double pschu = csch * real(s0);
+ double s0mag = cs0 * cabs(s0);
fprintf(
output,
" QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",
@@ -491,13 +496,13 @@ void sphere(string config_file, string data_file, string output_path) {
fprintf(output, " SPHERE %2d\n", ns);
fprintf(
output, " SAS(1,1)=%15.7lE%15.7lE, SAS(2,1)=%15.7lE%15.7lE\n",
- c1->sas[ns226][0][0].real(), c1->sas[ns226][0][0].imag(),
- c1->sas[ns226][1][0].real(), c1->sas[ns226][1][0].imag()
+ real(c1->sas[ns226][0][0]), imag(c1->sas[ns226][0][0]),
+ real(c1->sas[ns226][1][0]), imag(c1->sas[ns226][1][0])
);
fprintf(
output, " SAS(1,2)=%15.7lE%15.7lE, SAS(2,2)=%15.7lE%15.7lE\n",
- c1->sas[ns226][0][1].real(), c1->sas[ns226][0][1].imag(),
- c1->sas[ns226][1][1].real(), c1->sas[ns226][1][1].imag()
+ real(c1->sas[ns226][0][1]), imag(c1->sas[ns226][0][1]),
+ real(c1->sas[ns226][1][1]), imag(c1->sas[ns226][1][1])
);
if (jths == 1 && jphs == 1)
fprintf(
@@ -523,9 +528,9 @@ void sphere(string config_file, string data_file, string output_path) {
else fprintf(output, "\n");
}
for (int vi = 0; vi < 16; vi++) {
- double value = vint[vi].real();
+ double value = real(vint[vi]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = vint[vi].imag();
+ value = imag(vint[vi]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
for (int imul = 0; imul < 4; imul++) {
diff --git a/src/testing/test_TEDF.cpp b/src/testing/test_TEDF.cpp
index 98a681d3895a63659978cd130352083a7028c6e7..b790afd907df1dbb229e0f0ac4b214fd42ea7d5b 100644
--- a/src/testing/test_TEDF.cpp
+++ b/src/testing/test_TEDF.cpp
@@ -1,11 +1,16 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
//! \file test_TEDF.cpp
-#include <complex>
#include <cstdio>
#include <exception>
#include <hdf5.h>
#include <string>
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
+
#ifndef INCLUDE_ERRORS_H_
#include "../include/errors.h"
#endif
diff --git a/src/testing/test_TTMS.cpp b/src/testing/test_TTMS.cpp
index 7c55d16a7939a67ba211f90644b23ba2d329852e..77b8d7c22b1fdc59a73a7165ae2369c2e3c368ff 100644
--- a/src/testing/test_TTMS.cpp
+++ b/src/testing/test_TTMS.cpp
@@ -1,11 +1,16 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
//! \file test_TTMS.cpp
-#include <complex>
#include <cstdio>
#include <exception>
#include <hdf5.h>
#include <string>
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
+
#ifndef INCLUDE_ERRORS_H_
#include "../include/errors.h"
#endif
diff --git a/src/trapping/cfrfme.cpp b/src/trapping/cfrfme.cpp
index c4866e1f65073bcc2f61ec54d7a23a125ac315ca..2cd6ab9f33a227aff22d4255d4cdf48d01ccdb08 100644
--- a/src/trapping/cfrfme.cpp
+++ b/src/trapping/cfrfme.cpp
@@ -1,14 +1,19 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file cfrfme.cpp
*
* C++ implementation of FRFME functions.
*/
-#include <complex>
#include <cstdio>
#include <exception>
#include <fstream>
#include <regex>
#include <string>
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
+
#ifndef INCLUDE_PARSERS_H_
#include "../include/Parsers.h"
#endif
@@ -47,10 +52,10 @@ void frfme(string data_file, string output_path) {
Swap2 *tt2 = NULL;
char namef[7];
char more;
- complex<double> **w = NULL;
- complex<double> *wk = NULL;
- const complex<double> cc0(0.0, 0.0);
- const complex<double> uim(0.0, 1.0);
+ dcomplex **w = NULL;
+ dcomplex *wk = NULL;
+ const dcomplex cc0 = 0.0 + 0.0 * I;
+ const dcomplex uim = 0.0 + 1.0 * I;
int line_count = 0, last_read_line = 0;
regex re = regex("-?[0-9]+");
string *file_lines = load_file(data_file, &line_count);
@@ -255,7 +260,7 @@ void frfme(string data_file, string output_path) {
double *_xv = tfrfme->get_x();
double *_yv = tfrfme->get_y();
double *_zv = tfrfme->get_z();
- complex<double> **_wsum = tfrfme->get_matrix();
+ dcomplex **_wsum = tfrfme->get_matrix();
for (int i24 = nxshpo; i24 <= nxs; i24++) {
_xv[i24] = _xv[i24 - 1] + delxyz;
_xv[nxv - i24 - 1] = -_xv[i24];
@@ -307,15 +312,15 @@ void frfme(string data_file, string output_path) {
tt2->set_param("nrvc", 1.0 * nrvc);
tt2->write_binary(temp_name2, "HDF5");
for (int j80 = jlmf; j80 <= jlml; j80++) {
- complex<double> *tt1_wk = tt1->get_vector();
+ dcomplex *tt1_wk = tt1->get_vector();
int wk_index = 0;
// w matrix
if (w != NULL) {
for (int wi = nkv - 1; wi > -1; wi--) delete[] w[wi];
delete[] w;
}
- w = new complex<double>*[nkv];
- for (int wi = 0; wi < nkv; wi++) w[wi] = new complex<double>[nkv]();
+ w = new dcomplex*[nkv];
+ for (int wi = 0; wi < nkv; wi++) w[wi] = new dcomplex[nkv]();
for (int jy50 = 0; jy50 < nkv; jy50++) {
for (int jx50 = 0; jx50 < nkv; jx50++) {
for (int wi = 0; wi < nlmmt; wi++) wk[wi] = tt1_wk[wk_index++];
@@ -331,19 +336,19 @@ void frfme(string data_file, string output_path) {
for (int ix65 = 0; ix65 < nxv; ix65++) {
double x = _xv[ix65];
ixyz++;
- complex<double> sumy = cc0;
+ dcomplex sumy = cc0;
for (int jy60 = 0; jy60 < nkv; jy60++) {
double vky = vkv[jy60];
double vkx = vkv[nkv - 1];
double vkzf = vkzm[0][jy60];
- complex<double> phasf = exp(uim * (-vkx * x + vky * y + vkzf * z));
+ dcomplex phasf = cexp(uim * (-vkx * x + vky * y + vkzf * z));
double vkzl = vkzm[nkv - 1][jy60];
- complex<double> phasl = exp(uim * (vkx * x + vky * y + vkzl * z));
- complex<double> sumx = 0.5 * (w[0][jy60] * phasf + w[nkv - 1][jy60] * phasl);
+ dcomplex phasl = cexp(uim * (vkx * x + vky * y + vkzl * z));
+ dcomplex sumx = 0.5 * (w[0][jy60] * phasf + w[nkv - 1][jy60] * phasl);
for (int jx55 = 2; jx55 <= nks; jx55++) {
vkx = vkv[jx55 - 1];
double vkz = vkzm[jx55 - 1][jy60];
- complex<double> phas = exp(uim * (vkx * x + vky * y + vkz * z));
+ dcomplex phas = cexp(uim * (vkx * x + vky * y + vkz * z));
sumx += (w[jx55 - 1][jy60] * phas);
} // jx55 loop
if (jy60 == 0 || jy60 == nkv - 1) sumx *= 0.5;
diff --git a/src/trapping/clffft.cpp b/src/trapping/clffft.cpp
index 7e3619a252ada2d7102a4d88b3b2a2ba97e60f19..a8033d8098c08694464fe1ab1ac206c84fd18f04 100644
--- a/src/trapping/clffft.cpp
+++ b/src/trapping/clffft.cpp
@@ -1,14 +1,19 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file clffft.cpp
*
* \brief C++ implementation of LFFFT functions.
*/
-#include <complex>
#include <cstdio>
#include <exception>
#include <fstream>
#include <regex>
#include <string>
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
+
#ifndef INCLUDE_PARSERS_H_
#include "../include/Parsers.h"
#endif
@@ -37,18 +42,18 @@ using namespace std;
* \param output_path: `string` Directory to write the output files in.
*/
void lffft(string data_file, string output_path) {
- const complex<double> uim(0.0, 1.0);
+ const dcomplex uim = 0.0 + 1.0 * I;
const double sq2i = 1.0 / sqrt(2.0);
- const complex<double> sq2iti = sq2i * uim;
+ const dcomplex sq2iti = sq2i * uim;
fstream tlfff, tlfft;
double ****zpv = NULL;
- complex<double> *ac = NULL, *ws = NULL, *wsl = NULL;
- complex<double> **am0m = NULL;
- complex<double> **amd = NULL;
+ dcomplex *ac = NULL, *ws = NULL, *wsl = NULL;
+ dcomplex **am0m = NULL;
+ dcomplex **amd = NULL;
int **indam = NULL;
- complex<double> *tmsm = NULL, *tmse = NULL, **tms = NULL;
- complex<double> **_wsum = NULL;
+ dcomplex *tmsm = NULL, *tmse = NULL, **tms = NULL;
+ dcomplex **_wsum = NULL;
int jft, jss, jtw;
int is, le, nvam = 0;
double vks, exris;
@@ -79,55 +84,55 @@ void lffft(string data_file, string output_path) {
cil->nlem = le * (le + 2);
cil->nlemt = cil->nlem + cil->nlem;
if (is >= 2222) { // label 120
- tms = new complex<double>*[le];
- for (int ti = 0; ti < le; ti++) tms[ti] = new complex<double>[3]();
+ tms = new dcomplex*[le];
+ for (int ti = 0; ti < le; ti++) tms[ti] = new dcomplex[3]();
// QUESTION|WARNING: original code uses LM without defining it. Where does it come from?
int lm = le;
for (int i = 0; i < lm; i++) {
double vreal, vimag;
ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));
ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
- tms[i][0] = complex<double>(vreal, vimag);
+ tms[i][0] = vreal + vimag * I;
ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));
ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
- tms[i][1] = complex<double>(vreal, vimag);
+ tms[i][1] = vreal + vimag * I;
ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));
ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
- tms[i][2] = complex<double>(vreal, vimag);
+ tms[i][2] = vreal + vimag * I;
} // i loop
} else if (is >= 1111) { // label 125
- tmsm = new complex<double>[le]();
- tmse = new complex<double>[le]();
+ tmsm = new dcomplex[le]();
+ tmse = new dcomplex[le]();
for (int i = 0; i < le; i++) {
double vreal, vimag;
ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));
ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
- tmsm[i] = complex<double>(vreal, vimag);
+ tmsm[i] = vreal + vimag * I;
ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));
ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
- tmse[i] = complex<double>(vreal, vimag);
+ tmse[i] = vreal + vimag * I;
} // i loop
} else if (is >= 0) { // label 135
- am0m = new complex<double>*[cil->nlemt];
- for (int ai = 0; ai < cil->nlemt; ai++) am0m[ai] = new complex<double>[cil->nlemt]();
+ am0m = new dcomplex*[cil->nlemt];
+ for (int ai = 0; ai < cil->nlemt; ai++) am0m[ai] = new dcomplex[cil->nlemt]();
for (int i = 0; i < cil->nlemt; i++) {
for (int j = 0; j < cil->nlemt; j++) {
double vreal, vimag;
ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));
ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
- am0m[i][j] = complex<double>(vreal, vimag);
+ am0m[i][j] = vreal + vimag * I;
} // j loop
} // i loop
} else if (is < 0) {
nvam = le * le + (le * (le + 1) * (le * 2 + 1)) / 3;
- amd = new complex<double>*[nvam];
- for (int ai = 0; ai < nvam; ai++) amd[ai] = new complex<double>[4]();
+ amd = new dcomplex*[nvam];
+ for (int ai = 0; ai < nvam; ai++) amd[ai] = new dcomplex[4]();
for (int i = 0; i < nvam; i++) {
for (int j = 0; j < 4; j++) {
double vreal, vimag;
ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));
ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
- amd[i][j] = complex<double>(vreal, vimag);
+ amd[i][j] = vreal + vimag * I;
} // j loop
} // i loop
indam = new int*[le];
@@ -265,8 +270,8 @@ void lffft(string data_file, string output_path) {
// label 160
const int nlmm = lm * (lm + 2);
const int nlmmt = nlmm + nlmm;
- ws = new complex<double>[nlmmt]();
- if (lm > le) wsl = new complex<double>[nlmmt]();
+ ws = new dcomplex[nlmmt]();
+ if (lm > le) wsl = new dcomplex[nlmmt]();
// FORTRAN writes two output formatted files without opening them
// explicitly. It is assumed thay can be opened here.
string out66_name = output_path + "/c_out66.txt";
@@ -282,12 +287,10 @@ void lffft(string data_file, string output_path) {
//binary_input.read(reinterpret_cast<char *>(&vimag), sizeof(double));
int row = i;
int col = (nyv * nxv * iz475) + (nxv * iy475) + ix475;
- complex<double> value = _wsum[row][col];
+ dcomplex value = _wsum[row][col];
if (lm <= le) {
- //ws[i] = complex<double>(vreal, vimag);
ws[i] = value;
} else { // label 170
- //wsl[i] = complex<double>(vreal, vimag);
wsl[i] = value;
for (int i175 = 0; i175 < cil->nlem; i175++) {
int ie = i175 + cil->nlem;
@@ -301,21 +304,21 @@ void lffft(string data_file, string output_path) {
if (is != 2222) {
if (is != 1111) {
if (is > 0) { // Goes to 305
- ac = new complex<double>[cil->nlemt]();
+ ac = new dcomplex[cil->nlemt]();
camp(ac, am0m, ws, cil);
// Goes to 445
} else if (is < 0) { // Goes to 405
- ac = new complex<double>[cil->nlemt]();
+ ac = new dcomplex[cil->nlemt]();
czamp(ac, amd, indam, ws, cil);
// Goes to 445
}
} else {
- ac = new complex<double>[cil->nlemt]();
+ ac = new dcomplex[cil->nlemt]();
samp(ac, tmsm, tmse, ws, cil);
// Goes to 445
}
} else {
- ac = new complex<double>[cil->nlemt]();
+ ac = new dcomplex[cil->nlemt]();
sampoa(ac, tms, ws, cil);
// Goes to 445
}
diff --git a/test_data/README.md b/test_data/README.md
index b075b9e1ef924bca544715635cc6537cd8112578..4db4e5af7e9ab05240bbd9a95fb0d6f0efb0a80c 100644
--- a/test_data/README.md
+++ b/test_data/README.md
@@ -1,3 +1,5 @@
+Distributed under the terms of *GNU GPL-3.0-or-later*
+
# Folder instructions
This directory collects the input files to test the code.