231 changed files with 15251 additions and 5934 deletions
-
7resources/3rdparty/eigen/.hg_archival.txt
-
1resources/3rdparty/eigen/.hgtags
-
3resources/3rdparty/eigen/.krazy
-
655resources/3rdparty/eigen/COPYING.LGPL
-
3resources/3rdparty/eigen/COPYING.README
-
40resources/3rdparty/eigen/Eigen/Core
-
2resources/3rdparty/eigen/Eigen/Eigenvalues
-
26resources/3rdparty/eigen/Eigen/MetisSupport
-
2resources/3rdparty/eigen/Eigen/OrderingMethods
-
17resources/3rdparty/eigen/Eigen/SparseLU
-
7resources/3rdparty/eigen/Eigen/src/Cholesky/LDLT.h
-
24resources/3rdparty/eigen/Eigen/src/CholmodSupport/CholmodSupport.h
-
24resources/3rdparty/eigen/Eigen/src/Core/Array.h
-
64resources/3rdparty/eigen/Eigen/src/Core/ArrayWrapper.h
-
755resources/3rdparty/eigen/Eigen/src/Core/AssignEvaluator.h
-
2resources/3rdparty/eigen/Eigen/src/Core/Assign_MKL.h
-
40resources/3rdparty/eigen/Eigen/src/Core/Block.h
-
1299resources/3rdparty/eigen/Eigen/src/Core/CoreEvaluators.h
-
18resources/3rdparty/eigen/Eigen/src/Core/CwiseBinaryOp.h
-
96resources/3rdparty/eigen/Eigen/src/Core/CwiseNullaryOp.h
-
8resources/3rdparty/eigen/Eigen/src/Core/CwiseUnaryOp.h
-
28resources/3rdparty/eigen/Eigen/src/Core/DenseBase.h
-
12resources/3rdparty/eigen/Eigen/src/Core/DenseCoeffsBase.h
-
83resources/3rdparty/eigen/Eigen/src/Core/DenseStorage.h
-
37resources/3rdparty/eigen/Eigen/src/Core/Diagonal.h
-
16resources/3rdparty/eigen/Eigen/src/Core/DiagonalMatrix.h
-
4resources/3rdparty/eigen/Eigen/src/Core/DiagonalProduct.h
-
4resources/3rdparty/eigen/Eigen/src/Core/Dot.h
-
64resources/3rdparty/eigen/Eigen/src/Core/Functors.h
-
18resources/3rdparty/eigen/Eigen/src/Core/Fuzzy.h
-
2resources/3rdparty/eigen/Eigen/src/Core/GeneralProduct.h
-
8resources/3rdparty/eigen/Eigen/src/Core/Map.h
-
50resources/3rdparty/eigen/Eigen/src/Core/MapBase.h
-
47resources/3rdparty/eigen/Eigen/src/Core/MathFunctions.h
-
28resources/3rdparty/eigen/Eigen/src/Core/MatrixBase.h
-
5resources/3rdparty/eigen/Eigen/src/Core/NoAlias.h
-
22resources/3rdparty/eigen/Eigen/src/Core/PermutationMatrix.h
-
131resources/3rdparty/eigen/Eigen/src/Core/PlainObjectBase.h
-
45resources/3rdparty/eigen/Eigen/src/Core/Product.h
-
12resources/3rdparty/eigen/Eigen/src/Core/ProductBase.h
-
411resources/3rdparty/eigen/Eigen/src/Core/ProductEvaluators.h
-
4resources/3rdparty/eigen/Eigen/src/Core/Random.h
-
254resources/3rdparty/eigen/Eigen/src/Core/Ref.h
-
28resources/3rdparty/eigen/Eigen/src/Core/Replicate.h
-
8resources/3rdparty/eigen/Eigen/src/Core/Select.h
-
1resources/3rdparty/eigen/Eigen/src/Core/StableNorm.h
-
34resources/3rdparty/eigen/Eigen/src/Core/Swap.h
-
28resources/3rdparty/eigen/Eigen/src/Core/Transpose.h
-
18resources/3rdparty/eigen/Eigen/src/Core/Transpositions.h
-
7resources/3rdparty/eigen/Eigen/src/Core/TriangularMatrix.h
-
24resources/3rdparty/eigen/Eigen/src/Core/VectorBlock.h
-
12resources/3rdparty/eigen/Eigen/src/Core/Visitor.h
-
39resources/3rdparty/eigen/Eigen/src/Core/arch/NEON/PacketMath.h
-
78resources/3rdparty/eigen/Eigen/src/Core/arch/SSE/MathFunctions.h
-
8resources/3rdparty/eigen/Eigen/src/Core/arch/SSE/PacketMath.h
-
2resources/3rdparty/eigen/Eigen/src/Core/products/GeneralBlockPanelKernel.h
-
24resources/3rdparty/eigen/Eigen/src/Core/products/GeneralMatrixVector.h
-
20resources/3rdparty/eigen/Eigen/src/Core/products/TriangularMatrixMatrix_MKL.h
-
8resources/3rdparty/eigen/Eigen/src/Core/products/TriangularMatrixVector_MKL.h
-
11resources/3rdparty/eigen/Eigen/src/Core/util/Constants.h
-
4resources/3rdparty/eigen/Eigen/src/Core/util/ForwardDeclarations.h
-
14resources/3rdparty/eigen/Eigen/src/Core/util/Macros.h
-
3resources/3rdparty/eigen/Eigen/src/Core/util/StaticAssert.h
-
25resources/3rdparty/eigen/Eigen/src/Core/util/XprHelper.h
-
2resources/3rdparty/eigen/Eigen/src/Eigen2Support/Geometry/AlignedBox.h
-
2resources/3rdparty/eigen/Eigen/src/Eigen2Support/Geometry/AngleAxis.h
-
2resources/3rdparty/eigen/Eigen/src/Eigen2Support/Geometry/Hyperplane.h
-
2resources/3rdparty/eigen/Eigen/src/Eigen2Support/Geometry/ParametrizedLine.h
-
2resources/3rdparty/eigen/Eigen/src/Eigen2Support/Geometry/Quaternion.h
-
2resources/3rdparty/eigen/Eigen/src/Eigen2Support/Geometry/Rotation2D.h
-
2resources/3rdparty/eigen/Eigen/src/Eigen2Support/Geometry/RotationBase.h
-
2resources/3rdparty/eigen/Eigen/src/Eigen2Support/Geometry/Scaling.h
-
2resources/3rdparty/eigen/Eigen/src/Eigen2Support/Geometry/Transform.h
-
2resources/3rdparty/eigen/Eigen/src/Eigen2Support/Geometry/Translation.h
-
2resources/3rdparty/eigen/Eigen/src/Eigen2Support/LeastSquares.h
-
2resources/3rdparty/eigen/Eigen/src/Eigen2Support/SVD.h
-
18resources/3rdparty/eigen/Eigen/src/Eigenvalues/ComplexEigenSolver.h
-
58resources/3rdparty/eigen/Eigen/src/Eigenvalues/ComplexSchur.h
-
2resources/3rdparty/eigen/Eigen/src/Eigenvalues/ComplexSchur_MKL.h
-
21resources/3rdparty/eigen/Eigen/src/Eigenvalues/EigenSolver.h
-
339resources/3rdparty/eigen/Eigen/src/Eigenvalues/GeneralizedEigenSolver.h
-
618resources/3rdparty/eigen/Eigen/src/Eigenvalues/RealQZ.h
-
50resources/3rdparty/eigen/Eigen/src/Eigenvalues/RealSchur.h
-
11resources/3rdparty/eigen/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
-
2resources/3rdparty/eigen/Eigen/src/Eigenvalues/SelfAdjointEigenSolver_MKL.h
-
4resources/3rdparty/eigen/Eigen/src/Geometry/AngleAxis.h
-
4resources/3rdparty/eigen/Eigen/src/Geometry/Hyperplane.h
-
4resources/3rdparty/eigen/Eigen/src/Geometry/ParametrizedLine.h
-
4resources/3rdparty/eigen/Eigen/src/Geometry/Quaternion.h
-
6resources/3rdparty/eigen/Eigen/src/Geometry/Rotation2D.h
-
6resources/3rdparty/eigen/Eigen/src/Geometry/Scaling.h
-
18resources/3rdparty/eigen/Eigen/src/Geometry/Transform.h
-
25resources/3rdparty/eigen/Eigen/src/Geometry/Umeyama.h
-
6resources/3rdparty/eigen/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h
-
101resources/3rdparty/eigen/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h
-
26resources/3rdparty/eigen/Eigen/src/Jacobi/Jacobi.h
-
6resources/3rdparty/eigen/Eigen/src/MetisSupport/CMakeLists.txt
-
138resources/3rdparty/eigen/Eigen/src/MetisSupport/MetisSupport.h
-
1849resources/3rdparty/eigen/Eigen/src/OrderingMethods/Eigen_Colamd.h
-
158resources/3rdparty/eigen/Eigen/src/OrderingMethods/Ordering.h
@ -1,4 +1,5 @@ |
|||
repo: 8a21fd850624c931e448cbcfb38168cb2717c790 |
|||
node: 43d9075b23ef596ddf396101956d06f446fc0765 |
|||
branch: 3.1 |
|||
tag: 3.1.1 |
|||
node: 5945cb388ded120eb6dd3a1dfd2766b8e83237a4 |
|||
branch: default |
|||
latesttag: 3.1.0-rc2 |
|||
latesttagdistance: 147 |
@ -1,3 +0,0 @@ |
|||
SKIP /disabled/ |
|||
SKIP /bench/ |
|||
SKIP /build/ |
@ -1,165 +1,502 @@ |
|||
GNU LESSER GENERAL PUBLIC LICENSE |
|||
Version 3, 29 June 2007 |
|||
Version 2.1, February 1999 |
|||
|
|||
Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/> |
|||
Copyright (C) 1991, 1999 Free Software Foundation, Inc. |
|||
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA |
|||
Everyone is permitted to copy and distribute verbatim copies |
|||
of this license document, but changing it is not allowed. |
|||
|
|||
|
|||
This version of the GNU Lesser General Public License incorporates |
|||
the terms and conditions of version 3 of the GNU General Public |
|||
License, supplemented by the additional permissions listed below. |
|||
|
|||
0. Additional Definitions. |
|||
|
|||
As used herein, "this License" refers to version 3 of the GNU Lesser |
|||
General Public License, and the "GNU GPL" refers to version 3 of the GNU |
|||
General Public License. |
|||
|
|||
"The Library" refers to a covered work governed by this License, |
|||
other than an Application or a Combined Work as defined below. |
|||
|
|||
An "Application" is any work that makes use of an interface provided |
|||
by the Library, but which is not otherwise based on the Library. |
|||
Defining a subclass of a class defined by the Library is deemed a mode |
|||
of using an interface provided by the Library. |
|||
|
|||
A "Combined Work" is a work produced by combining or linking an |
|||
Application with the Library. The particular version of the Library |
|||
with which the Combined Work was made is also called the "Linked |
|||
Version". |
|||
|
|||
The "Minimal Corresponding Source" for a Combined Work means the |
|||
Corresponding Source for the Combined Work, excluding any source code |
|||
for portions of the Combined Work that, considered in isolation, are |
|||
based on the Application, and not on the Linked Version. |
|||
|
|||
The "Corresponding Application Code" for a Combined Work means the |
|||
object code and/or source code for the Application, including any data |
|||
and utility programs needed for reproducing the Combined Work from the |
|||
Application, but excluding the System Libraries of the Combined Work. |
|||
|
|||
1. Exception to Section 3 of the GNU GPL. |
|||
|
|||
You may convey a covered work under sections 3 and 4 of this License |
|||
without being bound by section 3 of the GNU GPL. |
|||
|
|||
2. Conveying Modified Versions. |
|||
|
|||
If you modify a copy of the Library, and, in your modifications, a |
|||
facility refers to a function or data to be supplied by an Application |
|||
that uses the facility (other than as an argument passed when the |
|||
facility is invoked), then you may convey a copy of the modified |
|||
version: |
|||
|
|||
a) under this License, provided that you make a good faith effort to |
|||
ensure that, in the event an Application does not supply the |
|||
function or data, the facility still operates, and performs |
|||
whatever part of its purpose remains meaningful, or |
|||
|
|||
b) under the GNU GPL, with none of the additional permissions of |
|||
this License applicable to that copy. |
|||
|
|||
3. Object Code Incorporating Material from Library Header Files. |
|||
|
|||
The object code form of an Application may incorporate material from |
|||
a header file that is part of the Library. You may convey such object |
|||
code under terms of your choice, provided that, if the incorporated |
|||
material is not limited to numerical parameters, data structure |
|||
layouts and accessors, or small macros, inline functions and templates |
|||
(ten or fewer lines in length), you do both of the following: |
|||
|
|||
a) Give prominent notice with each copy of the object code that the |
|||
Library is used in it and that the Library and its use are |
|||
covered by this License. |
|||
|
|||
b) Accompany the object code with a copy of the GNU GPL and this license |
|||
document. |
|||
|
|||
4. Combined Works. |
|||
|
|||
You may convey a Combined Work under terms of your choice that, |
|||
taken together, effectively do not restrict modification of the |
|||
portions of the Library contained in the Combined Work and reverse |
|||
engineering for debugging such modifications, if you also do each of |
|||
the following: |
|||
|
|||
a) Give prominent notice with each copy of the Combined Work that |
|||
the Library is used in it and that the Library and its use are |
|||
covered by this License. |
|||
|
|||
b) Accompany the Combined Work with a copy of the GNU GPL and this license |
|||
document. |
|||
|
|||
c) For a Combined Work that displays copyright notices during |
|||
execution, include the copyright notice for the Library among |
|||
these notices, as well as a reference directing the user to the |
|||
copies of the GNU GPL and this license document. |
|||
|
|||
d) Do one of the following: |
|||
|
|||
0) Convey the Minimal Corresponding Source under the terms of this |
|||
License, and the Corresponding Application Code in a form |
|||
suitable for, and under terms that permit, the user to |
|||
recombine or relink the Application with a modified version of |
|||
the Linked Version to produce a modified Combined Work, in the |
|||
manner specified by section 6 of the GNU GPL for conveying |
|||
Corresponding Source. |
|||
|
|||
1) Use a suitable shared library mechanism for linking with the |
|||
Library. A suitable mechanism is one that (a) uses at run time |
|||
a copy of the Library already present on the user's computer |
|||
system, and (b) will operate properly with a modified version |
|||
of the Library that is interface-compatible with the Linked |
|||
Version. |
|||
|
|||
e) Provide Installation Information, but only if you would otherwise |
|||
be required to provide such information under section 6 of the |
|||
GNU GPL, and only to the extent that such information is |
|||
necessary to install and execute a modified version of the |
|||
Combined Work produced by recombining or relinking the |
|||
Application with a modified version of the Linked Version. (If |
|||
you use option 4d0, the Installation Information must accompany |
|||
the Minimal Corresponding Source and Corresponding Application |
|||
Code. If you use option 4d1, you must provide the Installation |
|||
Information in the manner specified by section 6 of the GNU GPL |
|||
for conveying Corresponding Source.) |
|||
|
|||
5. Combined Libraries. |
|||
|
|||
You may place library facilities that are a work based on the |
|||
Library side by side in a single library together with other library |
|||
facilities that are not Applications and are not covered by this |
|||
License, and convey such a combined library under terms of your |
|||
choice, if you do both of the following: |
|||
|
|||
a) Accompany the combined library with a copy of the same work based |
|||
on the Library, uncombined with any other library facilities, |
|||
conveyed under the terms of this License. |
|||
|
|||
b) Give prominent notice with the combined library that part of it |
|||
is a work based on the Library, and explaining where to find the |
|||
accompanying uncombined form of the same work. |
|||
|
|||
6. Revised Versions of the GNU Lesser General Public License. |
|||
|
|||
The Free Software Foundation may publish revised and/or new versions |
|||
of the GNU Lesser 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 |
|||
Library as you received it specifies that a certain numbered version |
|||
of the GNU Lesser General Public License "or any later version" |
|||
applies to it, you have the option of following the terms and |
|||
conditions either of that published version or of any later version |
|||
published by the Free Software Foundation. If the Library as you |
|||
received it does not specify a version number of the GNU Lesser |
|||
General Public License, you may choose any version of the GNU Lesser |
|||
General Public License ever published by the Free Software Foundation. |
|||
|
|||
If the Library as you received it specifies that a proxy can decide |
|||
whether future versions of the GNU Lesser General Public License shall |
|||
apply, that proxy's public statement of acceptance of any version is |
|||
permanent authorization for you to choose that version for the |
|||
[This is the first released version of the Lesser GPL. It also counts |
|||
as the successor of the GNU Library Public License, version 2, hence |
|||
the version number 2.1.] |
|||
|
|||
Preamble |
|||
|
|||
The licenses for most software are designed to take away your |
|||
freedom to share and change it. By contrast, the GNU General Public |
|||
Licenses are intended to guarantee your freedom to share and change |
|||
free software--to make sure the software is free for all its users. |
|||
|
|||
This license, the Lesser General Public License, applies to some |
|||
specially designated software packages--typically libraries--of the |
|||
Free Software Foundation and other authors who decide to use it. You |
|||
can use it too, but we suggest you first think carefully about whether |
|||
this license or the ordinary General Public License is the better |
|||
strategy to use in any particular case, based on the explanations below. |
|||
|
|||
When we speak of free software, we are referring to freedom of use, |
|||
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 this service if you wish); that you receive source code or can get |
|||
it if you want it; that you can change the software and use pieces of |
|||
it in new free programs; and that you are informed that you can do |
|||
these things. |
|||
|
|||
To protect your rights, we need to make restrictions that forbid |
|||
distributors to deny you these rights or to ask you to surrender these |
|||
rights. These restrictions translate to certain responsibilities for |
|||
you if you distribute copies of the library or if you modify it. |
|||
|
|||
For example, if you distribute copies of the library, whether gratis |
|||
or for a fee, you must give the recipients all the rights that we gave |
|||
you. You must make sure that they, too, receive or can get the source |
|||
code. If you link other code with the library, you must provide |
|||
complete object files to the recipients, so that they can relink them |
|||
with the library after making changes to the library and recompiling |
|||
it. And you must show them these terms so they know their rights. |
|||
|
|||
We protect your rights with a two-step method: (1) we copyright the |
|||
library, and (2) we offer you this license, which gives you legal |
|||
permission to copy, distribute and/or modify the library. |
|||
|
|||
To protect each distributor, we want to make it very clear that |
|||
there is no warranty for the free library. Also, if the library is |
|||
modified by someone else and passed on, the recipients should know |
|||
that what they have is not the original version, so that the original |
|||
author's reputation will not be affected by problems that might be |
|||
introduced by others. |
|||
|
|||
Finally, software patents pose a constant threat to the existence of |
|||
any free program. We wish to make sure that a company cannot |
|||
effectively restrict the users of a free program by obtaining a |
|||
restrictive license from a patent holder. Therefore, we insist that |
|||
any patent license obtained for a version of the library must be |
|||
consistent with the full freedom of use specified in this license. |
|||
|
|||
Most GNU software, including some libraries, is covered by the |
|||
ordinary GNU General Public License. This license, the GNU Lesser |
|||
General Public License, applies to certain designated libraries, and |
|||
is quite different from the ordinary General Public License. We use |
|||
this license for certain libraries in order to permit linking those |
|||
libraries into non-free programs. |
|||
|
|||
When a program is linked with a library, whether statically or using |
|||
a shared library, the combination of the two is legally speaking a |
|||
combined work, a derivative of the original library. The ordinary |
|||
General Public License therefore permits such linking only if the |
|||
entire combination fits its criteria of freedom. The Lesser General |
|||
Public License permits more lax criteria for linking other code with |
|||
the library. |
|||
|
|||
We call this license the "Lesser" General Public License because it |
|||
does Less to protect the user's freedom than the ordinary General |
|||
Public License. It also provides other free software developers Less |
|||
of an advantage over competing non-free programs. These disadvantages |
|||
are the reason we use the ordinary General Public License for many |
|||
libraries. However, the Lesser license provides advantages in certain |
|||
special circumstances. |
|||
|
|||
For example, on rare occasions, there may be a special need to |
|||
encourage the widest possible use of a certain library, so that it becomes |
|||
a de-facto standard. To achieve this, non-free programs must be |
|||
allowed to use the library. A more frequent case is that a free |
|||
library does the same job as widely used non-free libraries. In this |
|||
case, there is little to gain by limiting the free library to free |
|||
software only, so we use the Lesser General Public License. |
|||
|
|||
In other cases, permission to use a particular library in non-free |
|||
programs enables a greater number of people to use a large body of |
|||
free software. For example, permission to use the GNU C Library in |
|||
non-free programs enables many more people to use the whole GNU |
|||
operating system, as well as its variant, the GNU/Linux operating |
|||
system. |
|||
|
|||
Although the Lesser General Public License is Less protective of the |
|||
users' freedom, it does ensure that the user of a program that is |
|||
linked with the Library has the freedom and the wherewithal to run |
|||
that program using a modified version of the Library. |
|||
|
|||
The precise terms and conditions for copying, distribution and |
|||
modification follow. Pay close attention to the difference between a |
|||
"work based on the library" and a "work that uses the library". The |
|||
former contains code derived from the library, whereas the latter must |
|||
be combined with the library in order to run. |
|||
|
|||
GNU LESSER GENERAL PUBLIC LICENSE |
|||
TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION |
|||
|
|||
0. This License Agreement applies to any software library or other |
|||
program which contains a notice placed by the copyright holder or |
|||
other authorized party saying it may be distributed under the terms of |
|||
this Lesser General Public License (also called "this License"). |
|||
Each licensee is addressed as "you". |
|||
|
|||
A "library" means a collection of software functions and/or data |
|||
prepared so as to be conveniently linked with application programs |
|||
(which use some of those functions and data) to form executables. |
|||
|
|||
The "Library", below, refers to any such software library or work |
|||
which has been distributed under these terms. A "work based on the |
|||
Library" means either the Library or any derivative work under |
|||
copyright law: that is to say, a work containing the Library or a |
|||
portion of it, either verbatim or with modifications and/or translated |
|||
straightforwardly into another language. (Hereinafter, translation is |
|||
included without limitation in the term "modification".) |
|||
|
|||
"Source code" for a work means the preferred form of the work for |
|||
making modifications to it. For a library, complete source code means |
|||
all the source code for all modules it contains, plus any associated |
|||
interface definition files, plus the scripts used to control compilation |
|||
and installation of the library. |
|||
|
|||
Activities other than copying, distribution and modification are not |
|||
covered by this License; they are outside its scope. The act of |
|||
running a program using the Library is not restricted, and output from |
|||
such a program is covered only if its contents constitute a work based |
|||
on the Library (independent of the use of the Library in a tool for |
|||
writing it). Whether that is true depends on what the Library does |
|||
and what the program that uses the Library does. |
|||
|
|||
1. You may copy and distribute verbatim copies of the Library's |
|||
complete source code as you receive it, in any medium, provided that |
|||
you conspicuously and appropriately publish on each copy an |
|||
appropriate copyright notice and disclaimer of warranty; keep intact |
|||
all the notices that refer to this License and to the absence of any |
|||
warranty; and distribute a copy of this License along with the |
|||
Library. |
|||
|
|||
You may charge a fee for the physical act of transferring a copy, |
|||
and you may at your option offer warranty protection in exchange for a |
|||
fee. |
|||
|
|||
2. You may modify your copy or copies of the Library or any portion |
|||
of it, thus forming a work based on the Library, and copy and |
|||
distribute such modifications or work under the terms of Section 1 |
|||
above, provided that you also meet all of these conditions: |
|||
|
|||
a) The modified work must itself be a software library. |
|||
|
|||
b) You must cause the files modified to carry prominent notices |
|||
stating that you changed the files and the date of any change. |
|||
|
|||
c) You must cause the whole of the work to be licensed at no |
|||
charge to all third parties under the terms of this License. |
|||
|
|||
d) If a facility in the modified Library refers to a function or a |
|||
table of data to be supplied by an application program that uses |
|||
the facility, other than as an argument passed when the facility |
|||
is invoked, then you must make a good faith effort to ensure that, |
|||
in the event an application does not supply such function or |
|||
table, the facility still operates, and performs whatever part of |
|||
its purpose remains meaningful. |
|||
|
|||
(For example, a function in a library to compute square roots has |
|||
a purpose that is entirely well-defined independent of the |
|||
application. Therefore, Subsection 2d requires that any |
|||
application-supplied function or table used by this function must |
|||
be optional: if the application does not supply it, the square |
|||
root function must still compute square roots.) |
|||
|
|||
These requirements apply to the modified work as a whole. If |
|||
identifiable sections of that work are not derived from the Library, |
|||
and can be reasonably considered independent and separate works in |
|||
themselves, then this License, and its terms, do not apply to those |
|||
sections when you distribute them as separate works. But when you |
|||
distribute the same sections as part of a whole which is a work based |
|||
on the Library, the distribution of the whole must be on the terms of |
|||
this License, whose permissions for other licensees extend to the |
|||
entire whole, and thus to each and every part regardless of who wrote |
|||
it. |
|||
|
|||
Thus, it is not the intent of this section to claim rights or contest |
|||
your rights to work written entirely by you; rather, the intent is to |
|||
exercise the right to control the distribution of derivative or |
|||
collective works based on the Library. |
|||
|
|||
In addition, mere aggregation of another work not based on the Library |
|||
with the Library (or with a work based on the Library) on a volume of |
|||
a storage or distribution medium does not bring the other work under |
|||
the scope of this License. |
|||
|
|||
3. You may opt to apply the terms of the ordinary GNU General Public |
|||
License instead of this License to a given copy of the Library. To do |
|||
this, you must alter all the notices that refer to this License, so |
|||
that they refer to the ordinary GNU General Public License, version 2, |
|||
instead of to this License. (If a newer version than version 2 of the |
|||
ordinary GNU General Public License has appeared, then you can specify |
|||
that version instead if you wish.) Do not make any other change in |
|||
these notices. |
|||
|
|||
Once this change is made in a given copy, it is irreversible for |
|||
that copy, so the ordinary GNU General Public License applies to all |
|||
subsequent copies and derivative works made from that copy. |
|||
|
|||
This option is useful when you wish to copy part of the code of |
|||
the Library into a program that is not a library. |
|||
|
|||
4. You may copy and distribute the Library (or a portion or |
|||
derivative of it, under Section 2) in object code or executable form |
|||
under the terms of Sections 1 and 2 above provided that you accompany |
|||
it with the complete corresponding machine-readable source code, which |
|||
must be distributed under the terms of Sections 1 and 2 above on a |
|||
medium customarily used for software interchange. |
|||
|
|||
If distribution of object code is made by offering access to copy |
|||
from a designated place, then offering equivalent access to copy the |
|||
source code from the same place satisfies the requirement to |
|||
distribute the source code, even though third parties are not |
|||
compelled to copy the source along with the object code. |
|||
|
|||
5. A program that contains no derivative of any portion of the |
|||
Library, but is designed to work with the Library by being compiled or |
|||
linked with it, is called a "work that uses the Library". Such a |
|||
work, in isolation, is not a derivative work of the Library, and |
|||
therefore falls outside the scope of this License. |
|||
|
|||
However, linking a "work that uses the Library" with the Library |
|||
creates an executable that is a derivative of the Library (because it |
|||
contains portions of the Library), rather than a "work that uses the |
|||
library". The executable is therefore covered by this License. |
|||
Section 6 states terms for distribution of such executables. |
|||
|
|||
When a "work that uses the Library" uses material from a header file |
|||
that is part of the Library, the object code for the work may be a |
|||
derivative work of the Library even though the source code is not. |
|||
Whether this is true is especially significant if the work can be |
|||
linked without the Library, or if the work is itself a library. The |
|||
threshold for this to be true is not precisely defined by law. |
|||
|
|||
If such an object file uses only numerical parameters, data |
|||
structure layouts and accessors, and small macros and small inline |
|||
functions (ten lines or less in length), then the use of the object |
|||
file is unrestricted, regardless of whether it is legally a derivative |
|||
work. (Executables containing this object code plus portions of the |
|||
Library will still fall under Section 6.) |
|||
|
|||
Otherwise, if the work is a derivative of the Library, you may |
|||
distribute the object code for the work under the terms of Section 6. |
|||
Any executables containing that work also fall under Section 6, |
|||
whether or not they are linked directly with the Library itself. |
|||
|
|||
6. As an exception to the Sections above, you may also combine or |
|||
link a "work that uses the Library" with the Library to produce a |
|||
work containing portions of the Library, and distribute that work |
|||
under terms of your choice, provided that the terms permit |
|||
modification of the work for the customer's own use and reverse |
|||
engineering for debugging such modifications. |
|||
|
|||
You must give prominent notice with each copy of the work that the |
|||
Library is used in it and that the Library and its use are covered by |
|||
this License. You must supply a copy of this License. If the work |
|||
during execution displays copyright notices, you must include the |
|||
copyright notice for the Library among them, as well as a reference |
|||
directing the user to the copy of this License. Also, you must do one |
|||
of these things: |
|||
|
|||
a) Accompany the work with the complete corresponding |
|||
machine-readable source code for the Library including whatever |
|||
changes were used in the work (which must be distributed under |
|||
Sections 1 and 2 above); and, if the work is an executable linked |
|||
with the Library, with the complete machine-readable "work that |
|||
uses the Library", as object code and/or source code, so that the |
|||
user can modify the Library and then relink to produce a modified |
|||
executable containing the modified Library. (It is understood |
|||
that the user who changes the contents of definitions files in the |
|||
Library will not necessarily be able to recompile the application |
|||
to use the modified definitions.) |
|||
|
|||
b) Use a suitable shared library mechanism for linking with the |
|||
Library. A suitable mechanism is one that (1) uses at run time a |
|||
copy of the library already present on the user's computer system, |
|||
rather than copying library functions into the executable, and (2) |
|||
will operate properly with a modified version of the library, if |
|||
the user installs one, as long as the modified version is |
|||
interface-compatible with the version that the work was made with. |
|||
|
|||
c) Accompany the work with a written offer, valid for at |
|||
least three years, to give the same user the materials |
|||
specified in Subsection 6a, above, for a charge no more |
|||
than the cost of performing this distribution. |
|||
|
|||
d) If distribution of the work is made by offering access to copy |
|||
from a designated place, offer equivalent access to copy the above |
|||
specified materials from the same place. |
|||
|
|||
e) Verify that the user has already received a copy of these |
|||
materials or that you have already sent this user a copy. |
|||
|
|||
For an executable, the required form of the "work that uses the |
|||
Library" must include any data and utility programs needed for |
|||
reproducing the executable from it. However, as a special exception, |
|||
the materials to be distributed need not include anything that is |
|||
normally distributed (in either source or binary form) with the major |
|||
components (compiler, kernel, and so on) of the operating system on |
|||
which the executable runs, unless that component itself accompanies |
|||
the executable. |
|||
|
|||
It may happen that this requirement contradicts the license |
|||
restrictions of other proprietary libraries that do not normally |
|||
accompany the operating system. Such a contradiction means you cannot |
|||
use both them and the Library together in an executable that you |
|||
distribute. |
|||
|
|||
7. You may place library facilities that are a work based on the |
|||
Library side-by-side in a single library together with other library |
|||
facilities not covered by this License, and distribute such a combined |
|||
library, provided that the separate distribution of the work based on |
|||
the Library and of the other library facilities is otherwise |
|||
permitted, and provided that you do these two things: |
|||
|
|||
a) Accompany the combined library with a copy of the same work |
|||
based on the Library, uncombined with any other library |
|||
facilities. This must be distributed under the terms of the |
|||
Sections above. |
|||
|
|||
b) Give prominent notice with the combined library of the fact |
|||
that part of it is a work based on the Library, and explaining |
|||
where to find the accompanying uncombined form of the same work. |
|||
|
|||
8. You may not copy, modify, sublicense, link with, or distribute |
|||
the Library except as expressly provided under this License. Any |
|||
attempt otherwise to copy, modify, sublicense, link with, or |
|||
distribute the Library is void, and will automatically terminate your |
|||
rights under this License. However, parties who have received copies, |
|||
or rights, from you under this License will not have their licenses |
|||
terminated so long as such parties remain in full compliance. |
|||
|
|||
9. You are not required to accept this License, since you have not |
|||
signed it. However, nothing else grants you permission to modify or |
|||
distribute the Library or its derivative works. These actions are |
|||
prohibited by law if you do not accept this License. Therefore, by |
|||
modifying or distributing the Library (or any work based on the |
|||
Library), you indicate your acceptance of this License to do so, and |
|||
all its terms and conditions for copying, distributing or modifying |
|||
the Library or works based on it. |
|||
|
|||
10. Each time you redistribute the Library (or any work based on the |
|||
Library), the recipient automatically receives a license from the |
|||
original licensor to copy, distribute, link with or modify the Library |
|||
subject to these terms and conditions. You may not impose any further |
|||
restrictions on the recipients' exercise of the rights granted herein. |
|||
You are not responsible for enforcing compliance by third parties with |
|||
this License. |
|||
|
|||
11. If, as a consequence of a court judgment or allegation of patent |
|||
infringement or for any other reason (not limited to patent issues), |
|||
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 |
|||
distribute so as to satisfy simultaneously your obligations under this |
|||
License and any other pertinent obligations, then as a consequence you |
|||
may not distribute the Library at all. For example, if a patent |
|||
license would not permit royalty-free redistribution of the Library by |
|||
all those who receive copies directly or indirectly through you, then |
|||
the only way you could satisfy both it and this License would be to |
|||
refrain entirely from distribution of the Library. |
|||
|
|||
If any portion of this section is held invalid or unenforceable under any |
|||
particular circumstance, the balance of the section is intended to apply, |
|||
and the section as a whole is intended to apply in other circumstances. |
|||
|
|||
It is not the purpose of this section to induce you to infringe any |
|||
patents or other property right claims or to contest validity of any |
|||
such claims; this section has the sole purpose of protecting the |
|||
integrity of the free software distribution system which is |
|||
implemented by public license practices. Many people have made |
|||
generous contributions to the wide range of software distributed |
|||
through that system in reliance on consistent application of that |
|||
system; it is up to the author/donor to decide if he or she is willing |
|||
to distribute software through any other system and a licensee cannot |
|||
impose that choice. |
|||
|
|||
This section is intended to make thoroughly clear what is believed to |
|||
be a consequence of the rest of this License. |
|||
|
|||
12. If the distribution and/or use of the Library is restricted in |
|||
certain countries either by patents or by copyrighted interfaces, the |
|||
original copyright holder who places the Library under this License may add |
|||
an explicit geographical distribution limitation excluding those countries, |
|||
so that distribution is permitted only in or among countries not thus |
|||
excluded. In such case, this License incorporates the limitation as if |
|||
written in the body of this License. |
|||
|
|||
13. The Free Software Foundation may publish revised and/or new |
|||
versions of the Lesser 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 Library |
|||
specifies a version number of this License which applies to it and |
|||
"any later version", you have the option of following the terms and |
|||
conditions either of that version or of any later version published by |
|||
the Free Software Foundation. If the Library does not specify a |
|||
license version number, you may choose any version ever published by |
|||
the Free Software Foundation. |
|||
|
|||
14. If you wish to incorporate parts of the Library into other free |
|||
programs whose distribution conditions are incompatible with these, |
|||
write to the author to ask for permission. For software which is |
|||
copyrighted by the Free Software Foundation, write to the Free |
|||
Software Foundation; we sometimes make exceptions for this. Our |
|||
decision will be guided by the two goals of preserving the free status |
|||
of all derivatives of our free software and of promoting the sharing |
|||
and reuse of software generally. |
|||
|
|||
NO WARRANTY |
|||
|
|||
15. BECAUSE THE LIBRARY IS LICENSED FREE OF CHARGE, THERE IS NO |
|||
WARRANTY FOR THE LIBRARY, TO THE EXTENT PERMITTED BY APPLICABLE LAW. |
|||
EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR |
|||
OTHER PARTIES PROVIDE THE LIBRARY "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 |
|||
LIBRARY IS WITH YOU. SHOULD THE LIBRARY PROVE DEFECTIVE, YOU ASSUME |
|||
THE COST OF ALL NECESSARY SERVICING, REPAIR OR CORRECTION. |
|||
|
|||
16. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN |
|||
WRITING WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY |
|||
AND/OR REDISTRIBUTE THE LIBRARY 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 |
|||
LIBRARY (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 LIBRARY TO OPERATE WITH ANY OTHER SOFTWARE), EVEN IF |
|||
SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH |
|||
DAMAGES. |
|||
|
|||
END OF TERMS AND CONDITIONS |
|||
|
|||
How to Apply These Terms to Your New Libraries |
|||
|
|||
If you develop a new library, and you want it to be of the greatest |
|||
possible use to the public, we recommend making it free software that |
|||
everyone can redistribute and change. You can do so by permitting |
|||
redistribution under these terms (or, alternatively, under the terms of the |
|||
ordinary General Public License). |
|||
|
|||
To apply these terms, attach the following notices to the library. It is |
|||
safest to attach them to the start of each source file to most effectively |
|||
convey 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 library's name and a brief idea of what it does.> |
|||
Copyright (C) <year> <name of author> |
|||
|
|||
This library is free software; you can redistribute it and/or |
|||
modify it under the terms of the GNU Lesser General Public |
|||
License as published by the Free Software Foundation; either |
|||
version 2.1 of the License, or (at your option) any later version. |
|||
|
|||
This library 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 |
|||
Lesser General Public License for more details. |
|||
|
|||
You should have received a copy of the GNU Lesser General Public |
|||
License along with this library; if not, write to the Free Software |
|||
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA |
|||
|
|||
Also add information on how to contact you by electronic and paper mail. |
|||
|
|||
You should also get your employer (if you work as a programmer) or your |
|||
school, if any, to sign a "copyright disclaimer" for the library, if |
|||
necessary. Here is a sample; alter the names: |
|||
|
|||
Yoyodyne, Inc., hereby disclaims all copyright interest in the |
|||
library `Frob' (a library for tweaking knobs) written by James Random Hacker. |
|||
|
|||
<signature of Ty Coon>, 1 April 1990 |
|||
Ty Coon, President of Vice |
|||
|
|||
That's all there is to it! |
@ -0,0 +1,26 @@ |
|||
#ifndef EIGEN_METISSUPPORT_MODULE_H |
|||
#define EIGEN_METISSUPPORT_MODULE_H |
|||
|
|||
#include "SparseCore" |
|||
|
|||
#include "src/Core/util/DisableStupidWarnings.h" |
|||
|
|||
extern "C" { |
|||
#include <metis.h> |
|||
} |
|||
|
|||
|
|||
/** \ingroup Support_modules |
|||
* \defgroup MetisSupport_Module MetisSupport module |
|||
* |
|||
* \code |
|||
* #include <Eigen/MetisSupport> |
|||
* \endcode |
|||
*/ |
|||
|
|||
|
|||
#include "src/MetisSupport/MetisSupport.h" |
|||
|
|||
#include "src/Core/util/ReenableStupidWarnings.h" |
|||
|
|||
#endif // EIGEN_METISSUPPORT_MODULE_H |
@ -0,0 +1,17 @@ |
|||
#ifndef EIGEN_SPARSELU_MODULE_H |
|||
#define EIGEN_SPARSELU_MODULE_H |
|||
|
|||
#include "SparseCore" |
|||
|
|||
|
|||
/** \ingroup Sparse_modules |
|||
* \defgroup SparseLU_Module SparseLU module |
|||
* |
|||
*/ |
|||
|
|||
// Ordering interface |
|||
#include "OrderingMethods" |
|||
|
|||
#include "src/SparseLU/SparseLU.h" |
|||
|
|||
#endif // EIGEN_SPARSELU_MODULE_H |
@ -0,0 +1,755 @@ |
|||
// This file is part of Eigen, a lightweight C++ template library |
|||
// for linear algebra. |
|||
// |
|||
// Copyright (C) 2011 Benoit Jacob <jacob.benoit.1@gmail.com> |
|||
// Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr> |
|||
// Copyright (C) 2011-2012 Jitse Niesen <jitse@maths.leeds.ac.uk> |
|||
// |
|||
// This Source Code Form is subject to the terms of the Mozilla |
|||
// Public License v. 2.0. If a copy of the MPL was not distributed |
|||
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. |
|||
|
|||
#ifndef EIGEN_ASSIGN_EVALUATOR_H |
|||
#define EIGEN_ASSIGN_EVALUATOR_H |
|||
|
|||
namespace Eigen { |
|||
|
|||
// This implementation is based on Assign.h |
|||
|
|||
namespace internal { |
|||
|
|||
/*************************************************************************** |
|||
* Part 1 : the logic deciding a strategy for traversal and unrolling * |
|||
***************************************************************************/ |
|||
|
|||
// copy_using_evaluator_traits is based on assign_traits |
|||
|
|||
template <typename Derived, typename OtherDerived> |
|||
struct copy_using_evaluator_traits |
|||
{ |
|||
public: |
|||
enum { |
|||
DstIsAligned = Derived::Flags & AlignedBit, |
|||
DstHasDirectAccess = Derived::Flags & DirectAccessBit, |
|||
SrcIsAligned = OtherDerived::Flags & AlignedBit, |
|||
JointAlignment = bool(DstIsAligned) && bool(SrcIsAligned) ? Aligned : Unaligned, |
|||
SrcEvalBeforeAssign = (evaluator_traits<OtherDerived>::HasEvalTo == 1) |
|||
}; |
|||
|
|||
private: |
|||
enum { |
|||
InnerSize = int(Derived::IsVectorAtCompileTime) ? int(Derived::SizeAtCompileTime) |
|||
: int(Derived::Flags)&RowMajorBit ? int(Derived::ColsAtCompileTime) |
|||
: int(Derived::RowsAtCompileTime), |
|||
InnerMaxSize = int(Derived::IsVectorAtCompileTime) ? int(Derived::MaxSizeAtCompileTime) |
|||
: int(Derived::Flags)&RowMajorBit ? int(Derived::MaxColsAtCompileTime) |
|||
: int(Derived::MaxRowsAtCompileTime), |
|||
MaxSizeAtCompileTime = Derived::SizeAtCompileTime, |
|||
PacketSize = packet_traits<typename Derived::Scalar>::size |
|||
}; |
|||
|
|||
enum { |
|||
StorageOrdersAgree = (int(Derived::IsRowMajor) == int(OtherDerived::IsRowMajor)), |
|||
MightVectorize = StorageOrdersAgree |
|||
&& (int(Derived::Flags) & int(OtherDerived::Flags) & ActualPacketAccessBit), |
|||
MayInnerVectorize = MightVectorize && int(InnerSize)!=Dynamic && int(InnerSize)%int(PacketSize)==0 |
|||
&& int(DstIsAligned) && int(SrcIsAligned), |
|||
MayLinearize = StorageOrdersAgree && (int(Derived::Flags) & int(OtherDerived::Flags) & LinearAccessBit), |
|||
MayLinearVectorize = MightVectorize && MayLinearize && DstHasDirectAccess |
|||
&& (DstIsAligned || MaxSizeAtCompileTime == Dynamic), |
|||
/* If the destination isn't aligned, we have to do runtime checks and we don't unroll, |
|||
so it's only good for large enough sizes. */ |
|||
MaySliceVectorize = MightVectorize && DstHasDirectAccess |
|||
&& (int(InnerMaxSize)==Dynamic || int(InnerMaxSize)>=3*PacketSize) |
|||
/* slice vectorization can be slow, so we only want it if the slices are big, which is |
|||
indicated by InnerMaxSize rather than InnerSize, think of the case of a dynamic block |
|||
in a fixed-size matrix */ |
|||
}; |
|||
|
|||
public: |
|||
enum { |
|||
Traversal = int(SrcEvalBeforeAssign) ? int(AllAtOnceTraversal) |
|||
: int(MayInnerVectorize) ? int(InnerVectorizedTraversal) |
|||
: int(MayLinearVectorize) ? int(LinearVectorizedTraversal) |
|||
: int(MaySliceVectorize) ? int(SliceVectorizedTraversal) |
|||
: int(MayLinearize) ? int(LinearTraversal) |
|||
: int(DefaultTraversal), |
|||
Vectorized = int(Traversal) == InnerVectorizedTraversal |
|||
|| int(Traversal) == LinearVectorizedTraversal |
|||
|| int(Traversal) == SliceVectorizedTraversal |
|||
}; |
|||
|
|||
private: |
|||
enum { |
|||
UnrollingLimit = EIGEN_UNROLLING_LIMIT * (Vectorized ? int(PacketSize) : 1), |
|||
MayUnrollCompletely = int(Derived::SizeAtCompileTime) != Dynamic |
|||
&& int(OtherDerived::CoeffReadCost) != Dynamic |
|||
&& int(Derived::SizeAtCompileTime) * int(OtherDerived::CoeffReadCost) <= int(UnrollingLimit), |
|||
MayUnrollInner = int(InnerSize) != Dynamic |
|||
&& int(OtherDerived::CoeffReadCost) != Dynamic |
|||
&& int(InnerSize) * int(OtherDerived::CoeffReadCost) <= int(UnrollingLimit) |
|||
}; |
|||
|
|||
public: |
|||
enum { |
|||
Unrolling = (int(Traversal) == int(InnerVectorizedTraversal) || int(Traversal) == int(DefaultTraversal)) |
|||
? ( |
|||
int(MayUnrollCompletely) ? int(CompleteUnrolling) |
|||
: int(MayUnrollInner) ? int(InnerUnrolling) |
|||
: int(NoUnrolling) |
|||
) |
|||
: int(Traversal) == int(LinearVectorizedTraversal) |
|||
? ( bool(MayUnrollCompletely) && bool(DstIsAligned) ? int(CompleteUnrolling) |
|||
: int(NoUnrolling) ) |
|||
: int(Traversal) == int(LinearTraversal) |
|||
? ( bool(MayUnrollCompletely) ? int(CompleteUnrolling) |
|||
: int(NoUnrolling) ) |
|||
: int(NoUnrolling) |
|||
}; |
|||
|
|||
#ifdef EIGEN_DEBUG_ASSIGN |
|||
static void debug() |
|||
{ |
|||
EIGEN_DEBUG_VAR(DstIsAligned) |
|||
EIGEN_DEBUG_VAR(SrcIsAligned) |
|||
EIGEN_DEBUG_VAR(JointAlignment) |
|||
EIGEN_DEBUG_VAR(InnerSize) |
|||
EIGEN_DEBUG_VAR(InnerMaxSize) |
|||
EIGEN_DEBUG_VAR(PacketSize) |
|||
EIGEN_DEBUG_VAR(StorageOrdersAgree) |
|||
EIGEN_DEBUG_VAR(MightVectorize) |
|||
EIGEN_DEBUG_VAR(MayLinearize) |
|||
EIGEN_DEBUG_VAR(MayInnerVectorize) |
|||
EIGEN_DEBUG_VAR(MayLinearVectorize) |
|||
EIGEN_DEBUG_VAR(MaySliceVectorize) |
|||
EIGEN_DEBUG_VAR(Traversal) |
|||
EIGEN_DEBUG_VAR(UnrollingLimit) |
|||
EIGEN_DEBUG_VAR(MayUnrollCompletely) |
|||
EIGEN_DEBUG_VAR(MayUnrollInner) |
|||
EIGEN_DEBUG_VAR(Unrolling) |
|||
} |
|||
#endif |
|||
}; |
|||
|
|||
/*************************************************************************** |
|||
* Part 2 : meta-unrollers |
|||
***************************************************************************/ |
|||
|
|||
/************************ |
|||
*** Default traversal *** |
|||
************************/ |
|||
|
|||
template<typename DstEvaluatorType, typename SrcEvaluatorType, int Index, int Stop> |
|||
struct copy_using_evaluator_DefaultTraversal_CompleteUnrolling |
|||
{ |
|||
typedef typename DstEvaluatorType::XprType DstXprType; |
|||
|
|||
enum { |
|||
outer = Index / DstXprType::InnerSizeAtCompileTime, |
|||
inner = Index % DstXprType::InnerSizeAtCompileTime |
|||
}; |
|||
|
|||
EIGEN_STRONG_INLINE static void run(DstEvaluatorType &dstEvaluator, |
|||
SrcEvaluatorType &srcEvaluator) |
|||
{ |
|||
dstEvaluator.copyCoeffByOuterInner(outer, inner, srcEvaluator); |
|||
copy_using_evaluator_DefaultTraversal_CompleteUnrolling |
|||
<DstEvaluatorType, SrcEvaluatorType, Index+1, Stop> |
|||
::run(dstEvaluator, srcEvaluator); |
|||
} |
|||
}; |
|||
|
|||
template<typename DstEvaluatorType, typename SrcEvaluatorType, int Stop> |
|||
struct copy_using_evaluator_DefaultTraversal_CompleteUnrolling<DstEvaluatorType, SrcEvaluatorType, Stop, Stop> |
|||
{ |
|||
EIGEN_STRONG_INLINE static void run(DstEvaluatorType&, SrcEvaluatorType&) { } |
|||
}; |
|||
|
|||
template<typename DstEvaluatorType, typename SrcEvaluatorType, int Index, int Stop> |
|||
struct copy_using_evaluator_DefaultTraversal_InnerUnrolling |
|||
{ |
|||
EIGEN_STRONG_INLINE static void run(DstEvaluatorType &dstEvaluator, |
|||
SrcEvaluatorType &srcEvaluator, |
|||
int outer) |
|||
{ |
|||
dstEvaluator.copyCoeffByOuterInner(outer, Index, srcEvaluator); |
|||
copy_using_evaluator_DefaultTraversal_InnerUnrolling |
|||
<DstEvaluatorType, SrcEvaluatorType, Index+1, Stop> |
|||
::run(dstEvaluator, srcEvaluator, outer); |
|||
} |
|||
}; |
|||
|
|||
template<typename DstEvaluatorType, typename SrcEvaluatorType, int Stop> |
|||
struct copy_using_evaluator_DefaultTraversal_InnerUnrolling<DstEvaluatorType, SrcEvaluatorType, Stop, Stop> |
|||
{ |
|||
EIGEN_STRONG_INLINE static void run(DstEvaluatorType&, SrcEvaluatorType&, int) { } |
|||
}; |
|||
|
|||
/*********************** |
|||
*** Linear traversal *** |
|||
***********************/ |
|||
|
|||
template<typename DstEvaluatorType, typename SrcEvaluatorType, int Index, int Stop> |
|||
struct copy_using_evaluator_LinearTraversal_CompleteUnrolling |
|||
{ |
|||
EIGEN_STRONG_INLINE static void run(DstEvaluatorType &dstEvaluator, |
|||
SrcEvaluatorType &srcEvaluator) |
|||
{ |
|||
dstEvaluator.copyCoeff(Index, srcEvaluator); |
|||
copy_using_evaluator_LinearTraversal_CompleteUnrolling |
|||
<DstEvaluatorType, SrcEvaluatorType, Index+1, Stop> |
|||
::run(dstEvaluator, srcEvaluator); |
|||
} |
|||
}; |
|||
|
|||
template<typename DstEvaluatorType, typename SrcEvaluatorType, int Stop> |
|||
struct copy_using_evaluator_LinearTraversal_CompleteUnrolling<DstEvaluatorType, SrcEvaluatorType, Stop, Stop> |
|||
{ |
|||
EIGEN_STRONG_INLINE static void run(DstEvaluatorType&, SrcEvaluatorType&) { } |
|||
}; |
|||
|
|||
/************************** |
|||
*** Inner vectorization *** |
|||
**************************/ |
|||
|
|||
template<typename DstEvaluatorType, typename SrcEvaluatorType, int Index, int Stop> |
|||
struct copy_using_evaluator_innervec_CompleteUnrolling |
|||
{ |
|||
typedef typename DstEvaluatorType::XprType DstXprType; |
|||
typedef typename SrcEvaluatorType::XprType SrcXprType; |
|||
|
|||
enum { |
|||
outer = Index / DstXprType::InnerSizeAtCompileTime, |
|||
inner = Index % DstXprType::InnerSizeAtCompileTime, |
|||
JointAlignment = copy_using_evaluator_traits<DstXprType,SrcXprType>::JointAlignment |
|||
}; |
|||
|
|||
EIGEN_STRONG_INLINE static void run(DstEvaluatorType &dstEvaluator, |
|||
SrcEvaluatorType &srcEvaluator) |
|||
{ |
|||
dstEvaluator.template copyPacketByOuterInner<Aligned, JointAlignment>(outer, inner, srcEvaluator); |
|||
enum { NextIndex = Index + packet_traits<typename DstXprType::Scalar>::size }; |
|||
copy_using_evaluator_innervec_CompleteUnrolling |
|||
<DstEvaluatorType, SrcEvaluatorType, NextIndex, Stop> |
|||
::run(dstEvaluator, srcEvaluator); |
|||
} |
|||
}; |
|||
|
|||
template<typename DstEvaluatorType, typename SrcEvaluatorType, int Stop> |
|||
struct copy_using_evaluator_innervec_CompleteUnrolling<DstEvaluatorType, SrcEvaluatorType, Stop, Stop> |
|||
{ |
|||
EIGEN_STRONG_INLINE static void run(DstEvaluatorType&, SrcEvaluatorType&) { } |
|||
}; |
|||
|
|||
template<typename DstEvaluatorType, typename SrcEvaluatorType, int Index, int Stop> |
|||
struct copy_using_evaluator_innervec_InnerUnrolling |
|||
{ |
|||
EIGEN_STRONG_INLINE static void run(DstEvaluatorType &dstEvaluator, |
|||
SrcEvaluatorType &srcEvaluator, |
|||
int outer) |
|||
{ |
|||
dstEvaluator.template copyPacketByOuterInner<Aligned, Aligned>(outer, Index, srcEvaluator); |
|||
typedef typename DstEvaluatorType::XprType DstXprType; |
|||
enum { NextIndex = Index + packet_traits<typename DstXprType::Scalar>::size }; |
|||
copy_using_evaluator_innervec_InnerUnrolling |
|||
<DstEvaluatorType, SrcEvaluatorType, NextIndex, Stop> |
|||
::run(dstEvaluator, srcEvaluator, outer); |
|||
} |
|||
}; |
|||
|
|||
template<typename DstEvaluatorType, typename SrcEvaluatorType, int Stop> |
|||
struct copy_using_evaluator_innervec_InnerUnrolling<DstEvaluatorType, SrcEvaluatorType, Stop, Stop> |
|||
{ |
|||
EIGEN_STRONG_INLINE static void run(DstEvaluatorType&, SrcEvaluatorType&, int) { } |
|||
}; |
|||
|
|||
/*************************************************************************** |
|||
* Part 3 : implementation of all cases |
|||
***************************************************************************/ |
|||
|
|||
// copy_using_evaluator_impl is based on assign_impl |
|||
|
|||
template<typename DstXprType, typename SrcXprType, |
|||
int Traversal = copy_using_evaluator_traits<DstXprType, SrcXprType>::Traversal, |
|||
int Unrolling = copy_using_evaluator_traits<DstXprType, SrcXprType>::Unrolling> |
|||
struct copy_using_evaluator_impl; |
|||
|
|||
/************************ |
|||
*** Default traversal *** |
|||
************************/ |
|||
|
|||
template<typename DstXprType, typename SrcXprType> |
|||
struct copy_using_evaluator_impl<DstXprType, SrcXprType, DefaultTraversal, NoUnrolling> |
|||
{ |
|||
static void run(DstXprType& dst, const SrcXprType& src) |
|||
{ |
|||
typedef typename evaluator<DstXprType>::type DstEvaluatorType; |
|||
typedef typename evaluator<SrcXprType>::type SrcEvaluatorType; |
|||
typedef typename DstXprType::Index Index; |
|||
|
|||
DstEvaluatorType dstEvaluator(dst); |
|||
SrcEvaluatorType srcEvaluator(src); |
|||
|
|||
for(Index outer = 0; outer < dst.outerSize(); ++outer) { |
|||
for(Index inner = 0; inner < dst.innerSize(); ++inner) { |
|||
dstEvaluator.copyCoeffByOuterInner(outer, inner, srcEvaluator); |
|||
} |
|||
} |
|||
} |
|||
}; |
|||
|
|||
template<typename DstXprType, typename SrcXprType> |
|||
struct copy_using_evaluator_impl<DstXprType, SrcXprType, DefaultTraversal, CompleteUnrolling> |
|||
{ |
|||
EIGEN_STRONG_INLINE static void run(DstXprType &dst, const SrcXprType &src) |
|||
{ |
|||
typedef typename evaluator<DstXprType>::type DstEvaluatorType; |
|||
typedef typename evaluator<SrcXprType>::type SrcEvaluatorType; |
|||
|
|||
DstEvaluatorType dstEvaluator(dst); |
|||
SrcEvaluatorType srcEvaluator(src); |
|||
|
|||
copy_using_evaluator_DefaultTraversal_CompleteUnrolling |
|||
<DstEvaluatorType, SrcEvaluatorType, 0, DstXprType::SizeAtCompileTime> |
|||
::run(dstEvaluator, srcEvaluator); |
|||
} |
|||
}; |
|||
|
|||
template<typename DstXprType, typename SrcXprType> |
|||
struct copy_using_evaluator_impl<DstXprType, SrcXprType, DefaultTraversal, InnerUnrolling> |
|||
{ |
|||
typedef typename DstXprType::Index Index; |
|||
EIGEN_STRONG_INLINE static void run(DstXprType &dst, const SrcXprType &src) |
|||
{ |
|||
typedef typename evaluator<DstXprType>::type DstEvaluatorType; |
|||
typedef typename evaluator<SrcXprType>::type SrcEvaluatorType; |
|||
|
|||
DstEvaluatorType dstEvaluator(dst); |
|||
SrcEvaluatorType srcEvaluator(src); |
|||
|
|||
const Index outerSize = dst.outerSize(); |
|||
for(Index outer = 0; outer < outerSize; ++outer) |
|||
copy_using_evaluator_DefaultTraversal_InnerUnrolling |
|||
<DstEvaluatorType, SrcEvaluatorType, 0, DstXprType::InnerSizeAtCompileTime> |
|||
::run(dstEvaluator, srcEvaluator, outer); |
|||
} |
|||
}; |
|||
|
|||
/*************************** |
|||
*** Linear vectorization *** |
|||
***************************/ |
|||
|
|||
template <bool IsAligned = false> |
|||
struct unaligned_copy_using_evaluator_impl |
|||
{ |
|||
// if IsAligned = true, then do nothing |
|||
template <typename SrcEvaluatorType, typename DstEvaluatorType> |
|||
static EIGEN_STRONG_INLINE void run(const SrcEvaluatorType&, DstEvaluatorType&, |
|||
typename SrcEvaluatorType::Index, typename SrcEvaluatorType::Index) {} |
|||
}; |
|||
|
|||
template <> |
|||
struct unaligned_copy_using_evaluator_impl<false> |
|||
{ |
|||
// MSVC must not inline this functions. If it does, it fails to optimize the |
|||
// packet access path. |
|||
#ifdef _MSC_VER |
|||
template <typename DstEvaluatorType, typename SrcEvaluatorType> |
|||
static EIGEN_DONT_INLINE void run(DstEvaluatorType &dstEvaluator, |
|||
const SrcEvaluatorType &srcEvaluator, |
|||
typename DstEvaluatorType::Index start, |
|||
typename DstEvaluatorType::Index end) |
|||
#else |
|||
template <typename DstEvaluatorType, typename SrcEvaluatorType> |
|||
static EIGEN_STRONG_INLINE void run(DstEvaluatorType &dstEvaluator, |
|||
const SrcEvaluatorType &srcEvaluator, |
|||
typename DstEvaluatorType::Index start, |
|||
typename DstEvaluatorType::Index end) |
|||
#endif |
|||
{ |
|||
for (typename DstEvaluatorType::Index index = start; index < end; ++index) |
|||
dstEvaluator.copyCoeff(index, srcEvaluator); |
|||
} |
|||
}; |
|||
|
|||
template<typename DstXprType, typename SrcXprType> |
|||
struct copy_using_evaluator_impl<DstXprType, SrcXprType, LinearVectorizedTraversal, NoUnrolling> |
|||
{ |
|||
EIGEN_STRONG_INLINE static void run(DstXprType &dst, const SrcXprType &src) |
|||
{ |
|||
typedef typename evaluator<DstXprType>::type DstEvaluatorType; |
|||
typedef typename evaluator<SrcXprType>::type SrcEvaluatorType; |
|||
typedef typename DstXprType::Index Index; |
|||
|
|||
DstEvaluatorType dstEvaluator(dst); |
|||
SrcEvaluatorType srcEvaluator(src); |
|||
|
|||
const Index size = dst.size(); |
|||
typedef packet_traits<typename DstXprType::Scalar> PacketTraits; |
|||
enum { |
|||
packetSize = PacketTraits::size, |
|||
dstIsAligned = int(copy_using_evaluator_traits<DstXprType,SrcXprType>::DstIsAligned), |
|||
dstAlignment = PacketTraits::AlignedOnScalar ? Aligned : dstIsAligned, |
|||
srcAlignment = copy_using_evaluator_traits<DstXprType,SrcXprType>::JointAlignment |
|||
}; |
|||
const Index alignedStart = dstIsAligned ? 0 : first_aligned(&dstEvaluator.coeffRef(0), size); |
|||
const Index alignedEnd = alignedStart + ((size-alignedStart)/packetSize)*packetSize; |
|||
|
|||
unaligned_copy_using_evaluator_impl<dstIsAligned!=0>::run(dstEvaluator, srcEvaluator, 0, alignedStart); |
|||
|
|||
for(Index index = alignedStart; index < alignedEnd; index += packetSize) |
|||
{ |
|||
dstEvaluator.template copyPacket<dstAlignment, srcAlignment>(index, srcEvaluator); |
|||
} |
|||
|
|||
unaligned_copy_using_evaluator_impl<>::run(dstEvaluator, srcEvaluator, alignedEnd, size); |
|||
} |
|||
}; |
|||
|
|||
template<typename DstXprType, typename SrcXprType> |
|||
struct copy_using_evaluator_impl<DstXprType, SrcXprType, LinearVectorizedTraversal, CompleteUnrolling> |
|||
{ |
|||
typedef typename DstXprType::Index Index; |
|||
EIGEN_STRONG_INLINE static void run(DstXprType &dst, const SrcXprType &src) |
|||
{ |
|||
typedef typename evaluator<DstXprType>::type DstEvaluatorType; |
|||
typedef typename evaluator<SrcXprType>::type SrcEvaluatorType; |
|||
|
|||
DstEvaluatorType dstEvaluator(dst); |
|||
SrcEvaluatorType srcEvaluator(src); |
|||
|
|||
enum { size = DstXprType::SizeAtCompileTime, |
|||
packetSize = packet_traits<typename DstXprType::Scalar>::size, |
|||
alignedSize = (size/packetSize)*packetSize }; |
|||
|
|||
copy_using_evaluator_innervec_CompleteUnrolling |
|||
<DstEvaluatorType, SrcEvaluatorType, 0, alignedSize> |
|||
::run(dstEvaluator, srcEvaluator); |
|||
copy_using_evaluator_DefaultTraversal_CompleteUnrolling |
|||
<DstEvaluatorType, SrcEvaluatorType, alignedSize, size> |
|||
::run(dstEvaluator, srcEvaluator); |
|||
} |
|||
}; |
|||
|
|||
/************************** |
|||
*** Inner vectorization *** |
|||
**************************/ |
|||
|
|||
template<typename DstXprType, typename SrcXprType> |
|||
struct copy_using_evaluator_impl<DstXprType, SrcXprType, InnerVectorizedTraversal, NoUnrolling> |
|||
{ |
|||
inline static void run(DstXprType &dst, const SrcXprType &src) |
|||
{ |
|||
typedef typename evaluator<DstXprType>::type DstEvaluatorType; |
|||
typedef typename evaluator<SrcXprType>::type SrcEvaluatorType; |
|||
typedef typename DstXprType::Index Index; |
|||
|
|||
DstEvaluatorType dstEvaluator(dst); |
|||
SrcEvaluatorType srcEvaluator(src); |
|||
|
|||
const Index innerSize = dst.innerSize(); |
|||
const Index outerSize = dst.outerSize(); |
|||
const Index packetSize = packet_traits<typename DstXprType::Scalar>::size; |
|||
for(Index outer = 0; outer < outerSize; ++outer) |
|||
for(Index inner = 0; inner < innerSize; inner+=packetSize) { |
|||
dstEvaluator.template copyPacketByOuterInner<Aligned, Aligned>(outer, inner, srcEvaluator); |
|||
} |
|||
} |
|||
}; |
|||
|
|||
template<typename DstXprType, typename SrcXprType> |
|||
struct copy_using_evaluator_impl<DstXprType, SrcXprType, InnerVectorizedTraversal, CompleteUnrolling> |
|||
{ |
|||
EIGEN_STRONG_INLINE static void run(DstXprType &dst, const SrcXprType &src) |
|||
{ |
|||
typedef typename evaluator<DstXprType>::type DstEvaluatorType; |
|||
typedef typename evaluator<SrcXprType>::type SrcEvaluatorType; |
|||
|
|||
DstEvaluatorType dstEvaluator(dst); |
|||
SrcEvaluatorType srcEvaluator(src); |
|||
|
|||
copy_using_evaluator_innervec_CompleteUnrolling |
|||
<DstEvaluatorType, SrcEvaluatorType, 0, DstXprType::SizeAtCompileTime> |
|||
::run(dstEvaluator, srcEvaluator); |
|||
} |
|||
}; |
|||
|
|||
template<typename DstXprType, typename SrcXprType> |
|||
struct copy_using_evaluator_impl<DstXprType, SrcXprType, InnerVectorizedTraversal, InnerUnrolling> |
|||
{ |
|||
typedef typename DstXprType::Index Index; |
|||
EIGEN_STRONG_INLINE static void run(DstXprType &dst, const SrcXprType &src) |
|||
{ |
|||
typedef typename evaluator<DstXprType>::type DstEvaluatorType; |
|||
typedef typename evaluator<SrcXprType>::type SrcEvaluatorType; |
|||
|
|||
DstEvaluatorType dstEvaluator(dst); |
|||
SrcEvaluatorType srcEvaluator(src); |
|||
|
|||
const Index outerSize = dst.outerSize(); |
|||
for(Index outer = 0; outer < outerSize; ++outer) |
|||
copy_using_evaluator_innervec_InnerUnrolling |
|||
<DstEvaluatorType, SrcEvaluatorType, 0, DstXprType::InnerSizeAtCompileTime> |
|||
::run(dstEvaluator, srcEvaluator, outer); |
|||
} |
|||
}; |
|||
|
|||
/*********************** |
|||
*** Linear traversal *** |
|||
***********************/ |
|||
|
|||
template<typename DstXprType, typename SrcXprType> |
|||
struct copy_using_evaluator_impl<DstXprType, SrcXprType, LinearTraversal, NoUnrolling> |
|||
{ |
|||
inline static void run(DstXprType &dst, const SrcXprType &src) |
|||
{ |
|||
typedef typename evaluator<DstXprType>::type DstEvaluatorType; |
|||
typedef typename evaluator<SrcXprType>::type SrcEvaluatorType; |
|||
typedef typename DstXprType::Index Index; |
|||
|
|||
DstEvaluatorType dstEvaluator(dst); |
|||
SrcEvaluatorType srcEvaluator(src); |
|||
|
|||
const Index size = dst.size(); |
|||
for(Index i = 0; i < size; ++i) |
|||
dstEvaluator.copyCoeff(i, srcEvaluator); |
|||
} |
|||
}; |
|||
|
|||
template<typename DstXprType, typename SrcXprType> |
|||
struct copy_using_evaluator_impl<DstXprType, SrcXprType, LinearTraversal, CompleteUnrolling> |
|||
{ |
|||
EIGEN_STRONG_INLINE static void run(DstXprType &dst, const SrcXprType &src) |
|||
{ |
|||
typedef typename evaluator<DstXprType>::type DstEvaluatorType; |
|||
typedef typename evaluator<SrcXprType>::type SrcEvaluatorType; |
|||
|
|||
DstEvaluatorType dstEvaluator(dst); |
|||
SrcEvaluatorType srcEvaluator(src); |
|||
|
|||
copy_using_evaluator_LinearTraversal_CompleteUnrolling |
|||
<DstEvaluatorType, SrcEvaluatorType, 0, DstXprType::SizeAtCompileTime> |
|||
::run(dstEvaluator, srcEvaluator); |
|||
} |
|||
}; |
|||
|
|||
/************************** |
|||
*** Slice vectorization *** |
|||
***************************/ |
|||
|
|||
template<typename DstXprType, typename SrcXprType> |
|||
struct copy_using_evaluator_impl<DstXprType, SrcXprType, SliceVectorizedTraversal, NoUnrolling> |
|||
{ |
|||
inline static void run(DstXprType &dst, const SrcXprType &src) |
|||
{ |
|||
typedef typename evaluator<DstXprType>::type DstEvaluatorType; |
|||
typedef typename evaluator<SrcXprType>::type SrcEvaluatorType; |
|||
typedef typename DstXprType::Index Index; |
|||
|
|||
DstEvaluatorType dstEvaluator(dst); |
|||
SrcEvaluatorType srcEvaluator(src); |
|||
|
|||
typedef packet_traits<typename DstXprType::Scalar> PacketTraits; |
|||
enum { |
|||
packetSize = PacketTraits::size, |
|||
alignable = PacketTraits::AlignedOnScalar, |
|||
dstAlignment = alignable ? Aligned : int(copy_using_evaluator_traits<DstXprType,SrcXprType>::DstIsAligned) |
|||
}; |
|||
const Index packetAlignedMask = packetSize - 1; |
|||
const Index innerSize = dst.innerSize(); |
|||
const Index outerSize = dst.outerSize(); |
|||
const Index alignedStep = alignable ? (packetSize - dst.outerStride() % packetSize) & packetAlignedMask : 0; |
|||
Index alignedStart = ((!alignable) || copy_using_evaluator_traits<DstXprType,SrcXprType>::DstIsAligned) ? 0 |
|||
: first_aligned(&dstEvaluator.coeffRef(0,0), innerSize); |
|||
|
|||
for(Index outer = 0; outer < outerSize; ++outer) |
|||
{ |
|||
const Index alignedEnd = alignedStart + ((innerSize-alignedStart) & ~packetAlignedMask); |
|||
// do the non-vectorizable part of the assignment |
|||
for(Index inner = 0; inner<alignedStart ; ++inner) { |
|||
dstEvaluator.copyCoeffByOuterInner(outer, inner, srcEvaluator); |
|||
} |
|||
|
|||
// do the vectorizable part of the assignment |
|||
for(Index inner = alignedStart; inner<alignedEnd; inner+=packetSize) { |
|||
dstEvaluator.template copyPacketByOuterInner<dstAlignment, Unaligned>(outer, inner, srcEvaluator); |
|||
} |
|||
|
|||
// do the non-vectorizable part of the assignment |
|||
for(Index inner = alignedEnd; inner<innerSize ; ++inner) { |
|||
dstEvaluator.copyCoeffByOuterInner(outer, inner, srcEvaluator); |
|||
} |
|||
|
|||
alignedStart = std::min<Index>((alignedStart+alignedStep)%packetSize, innerSize); |
|||
} |
|||
} |
|||
}; |
|||
|
|||
/**************************** |
|||
*** All-at-once traversal *** |
|||
****************************/ |
|||
|
|||
template<typename DstXprType, typename SrcXprType> |
|||
struct copy_using_evaluator_impl<DstXprType, SrcXprType, AllAtOnceTraversal, NoUnrolling> |
|||
{ |
|||
inline static void run(DstXprType &dst, const SrcXprType &src) |
|||
{ |
|||
typedef typename evaluator<DstXprType>::type DstEvaluatorType; |
|||
typedef typename evaluator<SrcXprType>::type SrcEvaluatorType; |
|||
typedef typename DstXprType::Index Index; |
|||
|
|||
DstEvaluatorType dstEvaluator(dst); |
|||
SrcEvaluatorType srcEvaluator(src); |
|||
|
|||
// Evaluate rhs in temporary to prevent aliasing problems in a = a * a; |
|||
// TODO: Do not pass the xpr object to evalTo() |
|||
srcEvaluator.evalTo(dstEvaluator, dst); |
|||
} |
|||
}; |
|||
|
|||
/*************************************************************************** |
|||
* Part 4 : Entry points |
|||
***************************************************************************/ |
|||
|
|||
// Based on DenseBase::LazyAssign() |
|||
|
|||
template<typename DstXprType, template <typename> class StorageBase, typename SrcXprType> |
|||
EIGEN_STRONG_INLINE |
|||
const DstXprType& copy_using_evaluator(const NoAlias<DstXprType, StorageBase>& dst, |
|||
const EigenBase<SrcXprType>& src) |
|||
{ |
|||
return noalias_copy_using_evaluator(dst.expression(), src.derived()); |
|||
} |
|||
|
|||
template<typename XprType, int AssumeAliasing = evaluator_traits<XprType>::AssumeAliasing> |
|||
struct AddEvalIfAssumingAliasing; |
|||
|
|||
template<typename XprType> |
|||
struct AddEvalIfAssumingAliasing<XprType, 0> |
|||
{ |
|||
static const XprType& run(const XprType& xpr) |
|||
{ |
|||
return xpr; |
|||
} |
|||
}; |
|||
|
|||
template<typename XprType> |
|||
struct AddEvalIfAssumingAliasing<XprType, 1> |
|||
{ |
|||
static const EvalToTemp<XprType> run(const XprType& xpr) |
|||
{ |
|||
return EvalToTemp<XprType>(xpr); |
|||
} |
|||
}; |
|||
|
|||
template<typename DstXprType, typename SrcXprType> |
|||
EIGEN_STRONG_INLINE |
|||
const DstXprType& copy_using_evaluator(const EigenBase<DstXprType>& dst, const EigenBase<SrcXprType>& src) |
|||
{ |
|||
return noalias_copy_using_evaluator(dst.const_cast_derived(), |
|||
AddEvalIfAssumingAliasing<SrcXprType>::run(src.derived())); |
|||
} |
|||
|
|||
template<typename DstXprType, typename SrcXprType> |
|||
EIGEN_STRONG_INLINE |
|||
const DstXprType& noalias_copy_using_evaluator(const PlainObjectBase<DstXprType>& dst, const EigenBase<SrcXprType>& src) |
|||
{ |
|||
#ifdef EIGEN_DEBUG_ASSIGN |
|||
internal::copy_using_evaluator_traits<DstXprType, SrcXprType>::debug(); |
|||
#endif |
|||
#ifdef EIGEN_NO_AUTOMATIC_RESIZING |
|||
eigen_assert((dst.size()==0 || (IsVectorAtCompileTime ? (dst.size() == src.size()) |
|||
: (dst.rows() == src.rows() && dst.cols() == src.cols()))) |
|||
&& "Size mismatch. Automatic resizing is disabled because EIGEN_NO_AUTOMATIC_RESIZING is defined"); |
|||
#else |
|||
dst.const_cast_derived().resizeLike(src.derived()); |
|||
#endif |
|||
return copy_using_evaluator_without_resizing(dst.const_cast_derived(), src.derived()); |
|||
} |
|||
|
|||
template<typename DstXprType, typename SrcXprType> |
|||
EIGEN_STRONG_INLINE |
|||
const DstXprType& noalias_copy_using_evaluator(const EigenBase<DstXprType>& dst, const EigenBase<SrcXprType>& src) |
|||
{ |
|||
return copy_using_evaluator_without_resizing(dst.const_cast_derived(), src.derived()); |
|||
} |
|||
|
|||
template<typename DstXprType, typename SrcXprType> |
|||
const DstXprType& copy_using_evaluator_without_resizing(const DstXprType& dst, const SrcXprType& src) |
|||
{ |
|||
#ifdef EIGEN_DEBUG_ASSIGN |
|||
internal::copy_using_evaluator_traits<DstXprType, SrcXprType>::debug(); |
|||
#endif |
|||
eigen_assert(dst.rows() == src.rows() && dst.cols() == src.cols()); |
|||
copy_using_evaluator_impl<DstXprType, SrcXprType>::run(const_cast<DstXprType&>(dst), src); |
|||
return dst; |
|||
} |
|||
|
|||
// Based on DenseBase::swap() |
|||
// TODO: Chech whether we need to do something special for swapping two |
|||
// Arrays or Matrices. |
|||
|
|||
template<typename DstXprType, typename SrcXprType> |
|||
void swap_using_evaluator(const DstXprType& dst, const SrcXprType& src) |
|||
{ |
|||
copy_using_evaluator(SwapWrapper<DstXprType>(const_cast<DstXprType&>(dst)), src); |
|||
} |
|||
|
|||
// Based on MatrixBase::operator+= (in CwiseBinaryOp.h) |
|||
template<typename DstXprType, typename SrcXprType> |
|||
void add_assign_using_evaluator(const MatrixBase<DstXprType>& dst, const MatrixBase<SrcXprType>& src) |
|||
{ |
|||
typedef typename DstXprType::Scalar Scalar; |
|||
SelfCwiseBinaryOp<internal::scalar_sum_op<Scalar>, DstXprType, SrcXprType> tmp(dst.const_cast_derived()); |
|||
copy_using_evaluator(tmp, src.derived()); |
|||
} |
|||
|
|||
// Based on ArrayBase::operator+= |
|||
template<typename DstXprType, typename SrcXprType> |
|||
void add_assign_using_evaluator(const ArrayBase<DstXprType>& dst, const ArrayBase<SrcXprType>& src) |
|||
{ |
|||
typedef typename DstXprType::Scalar Scalar; |
|||
SelfCwiseBinaryOp<internal::scalar_sum_op<Scalar>, DstXprType, SrcXprType> tmp(dst.const_cast_derived()); |
|||
copy_using_evaluator(tmp, src.derived()); |
|||
} |
|||
|
|||
// TODO: Add add_assign_using_evaluator for EigenBase ? |
|||
|
|||
template<typename DstXprType, typename SrcXprType> |
|||
void subtract_assign_using_evaluator(const MatrixBase<DstXprType>& dst, const MatrixBase<SrcXprType>& src) |
|||
{ |
|||
typedef typename DstXprType::Scalar Scalar; |
|||
SelfCwiseBinaryOp<internal::scalar_difference_op<Scalar>, DstXprType, SrcXprType> tmp(dst.const_cast_derived()); |
|||
copy_using_evaluator(tmp, src.derived()); |
|||
} |
|||
|
|||
template<typename DstXprType, typename SrcXprType> |
|||
void subtract_assign_using_evaluator(const ArrayBase<DstXprType>& dst, const ArrayBase<SrcXprType>& src) |
|||
{ |
|||
typedef typename DstXprType::Scalar Scalar; |
|||
SelfCwiseBinaryOp<internal::scalar_difference_op<Scalar>, DstXprType, SrcXprType> tmp(dst.const_cast_derived()); |
|||
copy_using_evaluator(tmp, src.derived()); |
|||
} |
|||
|
|||
template<typename DstXprType, typename SrcXprType> |
|||
void multiply_assign_using_evaluator(const ArrayBase<DstXprType>& dst, const ArrayBase<SrcXprType>& src) |
|||
{ |
|||
typedef typename DstXprType::Scalar Scalar; |
|||
SelfCwiseBinaryOp<internal::scalar_product_op<Scalar>, DstXprType, SrcXprType> tmp(dst.const_cast_derived()); |
|||
copy_using_evaluator(tmp, src.derived()); |
|||
} |
|||
|
|||
template<typename DstXprType, typename SrcXprType> |
|||
void divide_assign_using_evaluator(const ArrayBase<DstXprType>& dst, const ArrayBase<SrcXprType>& src) |
|||
{ |
|||
typedef typename DstXprType::Scalar Scalar; |
|||
SelfCwiseBinaryOp<internal::scalar_quotient_op<Scalar>, DstXprType, SrcXprType> tmp(dst.const_cast_derived()); |
|||
copy_using_evaluator(tmp, src.derived()); |
|||
} |
|||
|
|||
|
|||
} // namespace internal |
|||
|
|||
} // end namespace Eigen |
|||
|
|||
#endif // EIGEN_ASSIGN_EVALUATOR_H |
1299
resources/3rdparty/eigen/Eigen/src/Core/CoreEvaluators.h
File diff suppressed because it is too large
View File
File diff suppressed because it is too large
View File
@ -0,0 +1,411 @@ |
|||
// This file is part of Eigen, a lightweight C++ template library |
|||
// for linear algebra. |
|||
// |
|||
// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com> |
|||
// Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr> |
|||
// Copyright (C) 2011 Jitse Niesen <jitse@maths.leeds.ac.uk> |
|||
// |
|||
// This Source Code Form is subject to the terms of the Mozilla |
|||
// Public License v. 2.0. If a copy of the MPL was not distributed |
|||
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. |
|||
|
|||
|
|||
#ifndef EIGEN_PRODUCTEVALUATORS_H |
|||
#define EIGEN_PRODUCTEVALUATORS_H |
|||
|
|||
namespace Eigen { |
|||
|
|||
namespace internal { |
|||
|
|||
// We can evaluate the product either all at once, like GeneralProduct and its evalTo() function, or |
|||
// traverse the matrix coefficient by coefficient, like CoeffBasedProduct. Use the existing logic |
|||
// in ProductReturnType to decide. |
|||
|
|||
template<typename XprType, typename ProductType> |
|||
struct product_evaluator_dispatcher; |
|||
|
|||
template<typename Lhs, typename Rhs> |
|||
struct evaluator_impl<Product<Lhs, Rhs> > |
|||
: product_evaluator_dispatcher<Product<Lhs, Rhs>, typename ProductReturnType<Lhs, Rhs>::Type> |
|||
{ |
|||
typedef Product<Lhs, Rhs> XprType; |
|||
typedef product_evaluator_dispatcher<XprType, typename ProductReturnType<Lhs, Rhs>::Type> Base; |
|||
|
|||
evaluator_impl(const XprType& xpr) : Base(xpr) |
|||
{ } |
|||
}; |
|||
|
|||
template<typename XprType, typename ProductType> |
|||
struct product_evaluator_traits_dispatcher; |
|||
|
|||
template<typename Lhs, typename Rhs> |
|||
struct evaluator_traits<Product<Lhs, Rhs> > |
|||
: product_evaluator_traits_dispatcher<Product<Lhs, Rhs>, typename ProductReturnType<Lhs, Rhs>::Type> |
|||
{ |
|||
static const int AssumeAliasing = 1; |
|||
}; |
|||
|
|||
// Case 1: Evaluate all at once |
|||
// |
|||
// We can view the GeneralProduct class as a part of the product evaluator. |
|||
// Four sub-cases: InnerProduct, OuterProduct, GemmProduct and GemvProduct. |
|||
// InnerProduct is special because GeneralProduct does not have an evalTo() method in this case. |
|||
|
|||
template<typename Lhs, typename Rhs> |
|||
struct product_evaluator_traits_dispatcher<Product<Lhs, Rhs>, GeneralProduct<Lhs, Rhs, InnerProduct> > |
|||
{ |
|||
static const int HasEvalTo = 0; |
|||
}; |
|||
|
|||
template<typename Lhs, typename Rhs> |
|||
struct product_evaluator_dispatcher<Product<Lhs, Rhs>, GeneralProduct<Lhs, Rhs, InnerProduct> > |
|||
: public evaluator<typename Product<Lhs, Rhs>::PlainObject>::type |
|||
{ |
|||
typedef Product<Lhs, Rhs> XprType; |
|||
typedef typename XprType::PlainObject PlainObject; |
|||
typedef typename evaluator<PlainObject>::type evaluator_base; |
|||
|
|||
// TODO: Computation is too early (?) |
|||
product_evaluator_dispatcher(const XprType& xpr) : evaluator_base(m_result) |
|||
{ |
|||
m_result.coeffRef(0,0) = (xpr.lhs().transpose().cwiseProduct(xpr.rhs())).sum(); |
|||
} |
|||
|
|||
protected: |
|||
PlainObject m_result; |
|||
}; |
|||
|
|||
// For the other three subcases, simply call the evalTo() method of GeneralProduct |
|||
// TODO: GeneralProduct should take evaluators, not expression objects. |
|||
|
|||
template<typename Lhs, typename Rhs, int ProductType> |
|||
struct product_evaluator_traits_dispatcher<Product<Lhs, Rhs>, GeneralProduct<Lhs, Rhs, ProductType> > |
|||
{ |
|||
static const int HasEvalTo = 1; |
|||
}; |
|||
|
|||
template<typename Lhs, typename Rhs, int ProductType> |
|||
struct product_evaluator_dispatcher<Product<Lhs, Rhs>, GeneralProduct<Lhs, Rhs, ProductType> > |
|||
{ |
|||
typedef Product<Lhs, Rhs> XprType; |
|||
typedef typename XprType::PlainObject PlainObject; |
|||
typedef typename evaluator<PlainObject>::type evaluator_base; |
|||
|
|||
product_evaluator_dispatcher(const XprType& xpr) : m_xpr(xpr) |
|||
{ } |
|||
|
|||
template<typename DstEvaluatorType, typename DstXprType> |
|||
void evalTo(DstEvaluatorType /* not used */, DstXprType& dst) |
|||
{ |
|||
dst.resize(m_xpr.rows(), m_xpr.cols()); |
|||
GeneralProduct<Lhs, Rhs, ProductType>(m_xpr.lhs(), m_xpr.rhs()).evalTo(dst); |
|||
} |
|||
|
|||
protected: |
|||
const XprType& m_xpr; |
|||
}; |
|||
|
|||
// Case 2: Evaluate coeff by coeff |
|||
// |
|||
// This is mostly taken from CoeffBasedProduct.h |
|||
// The main difference is that we add an extra argument to the etor_product_*_impl::run() function |
|||
// for the inner dimension of the product, because evaluator object do not know their size. |
|||
|
|||
template<int Traversal, int UnrollingIndex, typename Lhs, typename Rhs, typename RetScalar> |
|||
struct etor_product_coeff_impl; |
|||
|
|||
template<int StorageOrder, int UnrollingIndex, typename Lhs, typename Rhs, typename Packet, int LoadMode> |
|||
struct etor_product_packet_impl; |
|||
|
|||
template<typename Lhs, typename Rhs, typename LhsNested, typename RhsNested, int Flags> |
|||
struct product_evaluator_traits_dispatcher<Product<Lhs, Rhs>, CoeffBasedProduct<LhsNested, RhsNested, Flags> > |
|||
{ |
|||
static const int HasEvalTo = 0; |
|||
}; |
|||
|
|||
template<typename Lhs, typename Rhs, typename LhsNested, typename RhsNested, int Flags> |
|||
struct product_evaluator_dispatcher<Product<Lhs, Rhs>, CoeffBasedProduct<LhsNested, RhsNested, Flags> > |
|||
: evaluator_impl_base<Product<Lhs, Rhs> > |
|||
{ |
|||
typedef Product<Lhs, Rhs> XprType; |
|||
typedef CoeffBasedProduct<LhsNested, RhsNested, Flags> CoeffBasedProductType; |
|||
|
|||
product_evaluator_dispatcher(const XprType& xpr) |
|||
: m_lhsImpl(xpr.lhs()), |
|||
m_rhsImpl(xpr.rhs()), |
|||
m_innerDim(xpr.lhs().cols()) |
|||
{ } |
|||
|
|||
typedef typename XprType::Index Index; |
|||
typedef typename XprType::Scalar Scalar; |
|||
typedef typename XprType::CoeffReturnType CoeffReturnType; |
|||
typedef typename XprType::PacketScalar PacketScalar; |
|||
typedef typename XprType::PacketReturnType PacketReturnType; |
|||
|
|||
// Everything below here is taken from CoeffBasedProduct.h |
|||
|
|||
enum { |
|||
RowsAtCompileTime = traits<CoeffBasedProductType>::RowsAtCompileTime, |
|||
PacketSize = packet_traits<Scalar>::size, |
|||
InnerSize = traits<CoeffBasedProductType>::InnerSize, |
|||
CoeffReadCost = traits<CoeffBasedProductType>::CoeffReadCost, |
|||
Unroll = CoeffReadCost != Dynamic && CoeffReadCost <= EIGEN_UNROLLING_LIMIT, |
|||
CanVectorizeInner = traits<CoeffBasedProductType>::CanVectorizeInner |
|||
}; |
|||
|
|||
typedef typename evaluator<Lhs>::type LhsEtorType; |
|||
typedef typename evaluator<Rhs>::type RhsEtorType; |
|||
typedef etor_product_coeff_impl<CanVectorizeInner ? InnerVectorizedTraversal : DefaultTraversal, |
|||
Unroll ? InnerSize-1 : Dynamic, |
|||
LhsEtorType, RhsEtorType, Scalar> CoeffImpl; |
|||
|
|||
const CoeffReturnType coeff(Index row, Index col) const |
|||
{ |
|||
Scalar res; |
|||
CoeffImpl::run(row, col, m_lhsImpl, m_rhsImpl, m_innerDim, res); |
|||
return res; |
|||
} |
|||
|
|||
/* Allow index-based non-packet access. It is impossible though to allow index-based packed access, |
|||
* which is why we don't set the LinearAccessBit. |
|||
*/ |
|||
const CoeffReturnType coeff(Index index) const |
|||
{ |
|||
Scalar res; |
|||
const Index row = RowsAtCompileTime == 1 ? 0 : index; |
|||
const Index col = RowsAtCompileTime == 1 ? index : 0; |
|||
CoeffImpl::run(row, col, m_lhsImpl, m_rhsImpl, m_innerDim, res); |
|||
return res; |
|||
} |
|||
|
|||
template<int LoadMode> |
|||
const PacketReturnType packet(Index row, Index col) const |
|||
{ |
|||
PacketScalar res; |
|||
typedef etor_product_packet_impl<Flags&RowMajorBit ? RowMajor : ColMajor, |
|||
Unroll ? InnerSize-1 : Dynamic, |
|||
LhsEtorType, RhsEtorType, PacketScalar, LoadMode> PacketImpl; |
|||
PacketImpl::run(row, col, m_lhsImpl, m_rhsImpl, m_innerDim, res); |
|||
return res; |
|||
} |
|||
|
|||
protected: |
|||
typename evaluator<Lhs>::type m_lhsImpl; |
|||
typename evaluator<Rhs>::type m_rhsImpl; |
|||
|
|||
// TODO: Get rid of m_innerDim if known at compile time |
|||
Index m_innerDim; |
|||
}; |
|||
|
|||
/*************************************************************************** |
|||
* Normal product .coeff() implementation (with meta-unrolling) |
|||
***************************************************************************/ |
|||
|
|||
/************************************** |
|||
*** Scalar path - no vectorization *** |
|||
**************************************/ |
|||
|
|||
template<int UnrollingIndex, typename Lhs, typename Rhs, typename RetScalar> |
|||
struct etor_product_coeff_impl<DefaultTraversal, UnrollingIndex, Lhs, Rhs, RetScalar> |
|||
{ |
|||
typedef typename Lhs::Index Index; |
|||
EIGEN_STRONG_INLINE static void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index innerDim, RetScalar &res) |
|||
{ |
|||
etor_product_coeff_impl<DefaultTraversal, UnrollingIndex-1, Lhs, Rhs, RetScalar>::run(row, col, lhs, rhs, innerDim, res); |
|||
res += lhs.coeff(row, UnrollingIndex) * rhs.coeff(UnrollingIndex, col); |
|||
} |
|||
}; |
|||
|
|||
template<typename Lhs, typename Rhs, typename RetScalar> |
|||
struct etor_product_coeff_impl<DefaultTraversal, 0, Lhs, Rhs, RetScalar> |
|||
{ |
|||
typedef typename Lhs::Index Index; |
|||
EIGEN_STRONG_INLINE static void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index /*innerDim*/, RetScalar &res) |
|||
{ |
|||
res = lhs.coeff(row, 0) * rhs.coeff(0, col); |
|||
} |
|||
}; |
|||
|
|||
template<typename Lhs, typename Rhs, typename RetScalar> |
|||
struct etor_product_coeff_impl<DefaultTraversal, Dynamic, Lhs, Rhs, RetScalar> |
|||
{ |
|||
typedef typename Lhs::Index Index; |
|||
EIGEN_STRONG_INLINE static void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index innerDim, RetScalar& res) |
|||
{ |
|||
eigen_assert(innerDim>0 && "you are using a non initialized matrix"); |
|||
res = lhs.coeff(row, 0) * rhs.coeff(0, col); |
|||
for(Index i = 1; i < innerDim; ++i) |
|||
res += lhs.coeff(row, i) * rhs.coeff(i, col); |
|||
} |
|||
}; |
|||
|
|||
/******************************************* |
|||
*** Scalar path with inner vectorization *** |
|||
*******************************************/ |
|||
|
|||
template<int UnrollingIndex, typename Lhs, typename Rhs, typename Packet> |
|||
struct etor_product_coeff_vectorized_unroller |
|||
{ |
|||
typedef typename Lhs::Index Index; |
|||
enum { PacketSize = packet_traits<typename Lhs::Scalar>::size }; |
|||
EIGEN_STRONG_INLINE static void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index innerDim, typename Lhs::PacketScalar &pres) |
|||
{ |
|||
etor_product_coeff_vectorized_unroller<UnrollingIndex-PacketSize, Lhs, Rhs, Packet>::run(row, col, lhs, rhs, innerDim, pres); |
|||
pres = padd(pres, pmul( lhs.template packet<Aligned>(row, UnrollingIndex) , rhs.template packet<Aligned>(UnrollingIndex, col) )); |
|||
} |
|||
}; |
|||
|
|||
template<typename Lhs, typename Rhs, typename Packet> |
|||
struct etor_product_coeff_vectorized_unroller<0, Lhs, Rhs, Packet> |
|||
{ |
|||
typedef typename Lhs::Index Index; |
|||
EIGEN_STRONG_INLINE static void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index /*innerDim*/, typename Lhs::PacketScalar &pres) |
|||
{ |
|||
pres = pmul(lhs.template packet<Aligned>(row, 0) , rhs.template packet<Aligned>(0, col)); |
|||
} |
|||
}; |
|||
|
|||
template<int UnrollingIndex, typename Lhs, typename Rhs, typename RetScalar> |
|||
struct etor_product_coeff_impl<InnerVectorizedTraversal, UnrollingIndex, Lhs, Rhs, RetScalar> |
|||
{ |
|||
typedef typename Lhs::PacketScalar Packet; |
|||
typedef typename Lhs::Index Index; |
|||
enum { PacketSize = packet_traits<typename Lhs::Scalar>::size }; |
|||
EIGEN_STRONG_INLINE static void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index innerDim, RetScalar &res) |
|||
{ |
|||
Packet pres; |
|||
etor_product_coeff_vectorized_unroller<UnrollingIndex+1-PacketSize, Lhs, Rhs, Packet>::run(row, col, lhs, rhs, innerDim, pres); |
|||
etor_product_coeff_impl<DefaultTraversal,UnrollingIndex,Lhs,Rhs,RetScalar>::run(row, col, lhs, rhs, innerDim, res); |
|||
res = predux(pres); |
|||
} |
|||
}; |
|||
|
|||
template<typename Lhs, typename Rhs, int LhsRows = Lhs::RowsAtCompileTime, int RhsCols = Rhs::ColsAtCompileTime> |
|||
struct etor_product_coeff_vectorized_dyn_selector |
|||
{ |
|||
typedef typename Lhs::Index Index; |
|||
EIGEN_STRONG_INLINE static void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index /*innerDim*/, typename Lhs::Scalar &res) |
|||
{ |
|||
res = lhs.row(row).transpose().cwiseProduct(rhs.col(col)).sum(); |
|||
} |
|||
}; |
|||
|
|||
// NOTE the 3 following specializations are because taking .col(0) on a vector is a bit slower |
|||
// NOTE maybe they are now useless since we have a specialization for Block<Matrix> |
|||
template<typename Lhs, typename Rhs, int RhsCols> |
|||
struct etor_product_coeff_vectorized_dyn_selector<Lhs,Rhs,1,RhsCols> |
|||
{ |
|||
typedef typename Lhs::Index Index; |
|||
EIGEN_STRONG_INLINE static void run(Index /*row*/, Index col, const Lhs& lhs, const Rhs& rhs, Index /*innerDim*/, typename Lhs::Scalar &res) |
|||
{ |
|||
res = lhs.transpose().cwiseProduct(rhs.col(col)).sum(); |
|||
} |
|||
}; |
|||
|
|||
template<typename Lhs, typename Rhs, int LhsRows> |
|||
struct etor_product_coeff_vectorized_dyn_selector<Lhs,Rhs,LhsRows,1> |
|||
{ |
|||
typedef typename Lhs::Index Index; |
|||
EIGEN_STRONG_INLINE static void run(Index row, Index /*col*/, const Lhs& lhs, const Rhs& rhs, Index /*innerDim*/, typename Lhs::Scalar &res) |
|||
{ |
|||
res = lhs.row(row).transpose().cwiseProduct(rhs).sum(); |
|||
} |
|||
}; |
|||
|
|||
template<typename Lhs, typename Rhs> |
|||
struct etor_product_coeff_vectorized_dyn_selector<Lhs,Rhs,1,1> |
|||
{ |
|||
typedef typename Lhs::Index Index; |
|||
EIGEN_STRONG_INLINE static void run(Index /*row*/, Index /*col*/, const Lhs& lhs, const Rhs& rhs, Index /*innerDim*/, typename Lhs::Scalar &res) |
|||
{ |
|||
res = lhs.transpose().cwiseProduct(rhs).sum(); |
|||
} |
|||
}; |
|||
|
|||
template<typename Lhs, typename Rhs, typename RetScalar> |
|||
struct etor_product_coeff_impl<InnerVectorizedTraversal, Dynamic, Lhs, Rhs, RetScalar> |
|||
{ |
|||
typedef typename Lhs::Index Index; |
|||
EIGEN_STRONG_INLINE static void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index innerDim, typename Lhs::Scalar &res) |
|||
{ |
|||
etor_product_coeff_vectorized_dyn_selector<Lhs,Rhs>::run(row, col, lhs, rhs, innerDim, res); |
|||
} |
|||
}; |
|||
|
|||
/******************* |
|||
*** Packet path *** |
|||
*******************/ |
|||
|
|||
template<int UnrollingIndex, typename Lhs, typename Rhs, typename Packet, int LoadMode> |
|||
struct etor_product_packet_impl<RowMajor, UnrollingIndex, Lhs, Rhs, Packet, LoadMode> |
|||
{ |
|||
typedef typename Lhs::Index Index; |
|||
EIGEN_STRONG_INLINE static void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index innerDim, Packet &res) |
|||
{ |
|||
etor_product_packet_impl<RowMajor, UnrollingIndex-1, Lhs, Rhs, Packet, LoadMode>::run(row, col, lhs, rhs, innerDim, res); |
|||
res = pmadd(pset1<Packet>(lhs.coeff(row, UnrollingIndex)), rhs.template packet<LoadMode>(UnrollingIndex, col), res); |
|||
} |
|||
}; |
|||
|
|||
template<int UnrollingIndex, typename Lhs, typename Rhs, typename Packet, int LoadMode> |
|||
struct etor_product_packet_impl<ColMajor, UnrollingIndex, Lhs, Rhs, Packet, LoadMode> |
|||
{ |
|||
typedef typename Lhs::Index Index; |
|||
EIGEN_STRONG_INLINE static void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index innerDim, Packet &res) |
|||
{ |
|||
etor_product_packet_impl<ColMajor, UnrollingIndex-1, Lhs, Rhs, Packet, LoadMode>::run(row, col, lhs, rhs, innerDim, res); |
|||
res = pmadd(lhs.template packet<LoadMode>(row, UnrollingIndex), pset1<Packet>(rhs.coeff(UnrollingIndex, col)), res); |
|||
} |
|||
}; |
|||
|
|||
template<typename Lhs, typename Rhs, typename Packet, int LoadMode> |
|||
struct etor_product_packet_impl<RowMajor, 0, Lhs, Rhs, Packet, LoadMode> |
|||
{ |
|||
typedef typename Lhs::Index Index; |
|||
EIGEN_STRONG_INLINE static void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index /*innerDim*/, Packet &res) |
|||
{ |
|||
res = pmul(pset1<Packet>(lhs.coeff(row, 0)),rhs.template packet<LoadMode>(0, col)); |
|||
} |
|||
}; |
|||
|
|||
template<typename Lhs, typename Rhs, typename Packet, int LoadMode> |
|||
struct etor_product_packet_impl<ColMajor, 0, Lhs, Rhs, Packet, LoadMode> |
|||
{ |
|||
typedef typename Lhs::Index Index; |
|||
EIGEN_STRONG_INLINE static void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index /*innerDim*/, Packet &res) |
|||
{ |
|||
res = pmul(lhs.template packet<LoadMode>(row, 0), pset1<Packet>(rhs.coeff(0, col))); |
|||
} |
|||
}; |
|||
|
|||
template<typename Lhs, typename Rhs, typename Packet, int LoadMode> |
|||
struct etor_product_packet_impl<RowMajor, Dynamic, Lhs, Rhs, Packet, LoadMode> |
|||
{ |
|||
typedef typename Lhs::Index Index; |
|||
EIGEN_STRONG_INLINE static void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index innerDim, Packet& res) |
|||
{ |
|||
eigen_assert(innerDim>0 && "you are using a non initialized matrix"); |
|||
res = pmul(pset1<Packet>(lhs.coeff(row, 0)),rhs.template packet<LoadMode>(0, col)); |
|||
for(Index i = 1; i < innerDim; ++i) |
|||
res = pmadd(pset1<Packet>(lhs.coeff(row, i)), rhs.template packet<LoadMode>(i, col), res); |
|||
} |
|||
}; |
|||
|
|||
template<typename Lhs, typename Rhs, typename Packet, int LoadMode> |
|||
struct etor_product_packet_impl<ColMajor, Dynamic, Lhs, Rhs, Packet, LoadMode> |
|||
{ |
|||
typedef typename Lhs::Index Index; |
|||
EIGEN_STRONG_INLINE static void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index innerDim, Packet& res) |
|||
{ |
|||
eigen_assert(innerDim>0 && "you are using a non initialized matrix"); |
|||
res = pmul(lhs.template packet<LoadMode>(row, 0), pset1<Packet>(rhs.coeff(0, col))); |
|||
for(Index i = 1; i < innerDim; ++i) |
|||
res = pmadd(lhs.template packet<LoadMode>(row, i), pset1<Packet>(rhs.coeff(i, col)), res); |
|||
} |
|||
}; |
|||
|
|||
} // end namespace internal |
|||
|
|||
} // end namespace Eigen |
|||
|
|||
#endif // EIGEN_PRODUCT_EVALUATORS_H |
@ -0,0 +1,254 @@ |
|||
// This file is part of Eigen, a lightweight C++ template library |
|||
// for linear algebra. |
|||
// |
|||
// Copyright (C) 2012 Gael Guennebaud <gael.guennebaud@inria.fr> |
|||
// |
|||
// This Source Code Form is subject to the terms of the Mozilla |
|||
// Public License v. 2.0. If a copy of the MPL was not distributed |
|||
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. |
|||
|
|||
#ifndef EIGEN_REF_H |
|||
#define EIGEN_REF_H |
|||
|
|||
namespace Eigen { |
|||
|
|||
template<typename Derived> class RefBase; |
|||
template<typename PlainObjectType, int Options = 0, |
|||
typename StrideType = typename internal::conditional<PlainObjectType::IsVectorAtCompileTime,InnerStride<1>,OuterStride<> >::type > class Ref; |
|||
|
|||
/** \class Ref |
|||
* \ingroup Core_Module |
|||
* |
|||
* \brief A matrix or vector expression mapping an existing expressions |
|||
* |
|||
* \tparam PlainObjectType the equivalent matrix type of the mapped data |
|||
* \tparam Options specifies whether the pointer is \c #Aligned, or \c #Unaligned. |
|||
* The default is \c #Unaligned. |
|||
* \tparam StrideType optionally specifies strides. By default, Ref implies a contiguous storage along the inner dimension (inner stride==1), |
|||
* but accept a variable outer stride (leading dimension). |
|||
* This can be overridden by specifying strides. |
|||
* The type passed here must be a specialization of the Stride template, see examples below. |
|||
* |
|||
* This class permits to write non template functions taking Eigen's object as parameters while limiting the number of copies. |
|||
* A Ref<> object can represent either a const expression or a l-value: |
|||
* \code |
|||
* // in-out argument: |
|||
* void foo1(Ref<VectorXf> x); |
|||
* |
|||
* // read-only const argument: |
|||
* void foo2(const Ref<const VectorXf>& x); |
|||
* \endcode |
|||
* |
|||
* In the in-out case, the input argument must satisfies the constraints of the actual Ref<> type, otherwise a compilation issue will be triggered. |
|||
* By default, a Ref<VectorXf> can reference any dense vector expression of float having a contiguous memory layout. |
|||
* Likewise, a Ref<MatrixXf> can reference any column major dense matrix expression of float whose column's elements are contiguously stored with |
|||
* the possibility to have a constant space inbetween each column, i.e.: the inner stride mmust be equal to 1, but the outer-stride (or leading dimension), |
|||
* can be greater than the number of rows. |
|||
* |
|||
* In the const case, if the input expression does not match the above requirement, then it is evaluated into a temporary before being passed to the function. |
|||
* Here are some examples: |
|||
* \code |
|||
* MatrixXf A; |
|||
* VectorXf a; |
|||
* foo1(a.head()); // OK |
|||
* foo1(A.col()); // OK |
|||
* foo1(A.row()); // compilation error because here innerstride!=1 |
|||
* foo2(A.row()); // The row is copied into a contiguous temporary |
|||
* foo2(2*a); // The expression is evaluated into a temporary |
|||
* foo2(A.col().segment(2,4)); // No temporary |
|||
* \endcode |
|||
* |
|||
* The range of inputs that can be referenced without temporary can be enlarged using the last two template parameter. |
|||
* Here is an example accepting an innerstride!=1: |
|||
* \code |
|||
* // in-out argument: |
|||
* void foo3(Ref<VectorXf,0,InnerStride<> > x); |
|||
* foo3(A.row()); // OK |
|||
* \endcode |
|||
* The downside here is that the function foo3 might be significantly slower than foo1 because it won't be able to exploit vectorization, and will involved more |
|||
* expensive address computations even if the input is contiguously stored in memory. To overcome this issue, one might propose to overloads internally calling a |
|||
* template function, e.g.: |
|||
* \code |
|||
* // in the .h: |
|||
* void foo(const Ref<MatrixXf>& A); |
|||
* void foo(const Ref<MatrixXf,0,Stride<> >& A); |
|||
* |
|||
* // in the .cpp: |
|||
* template<typename TypeOfA> void foo_impl(const TypeOfA& A) { |
|||
* ... // crazy code goes here |
|||
* } |
|||
* void foo(const Ref<MatrixXf>& A) { foo_impl(A); } |
|||
* void foo(const Ref<MatrixXf,0,Stride<> >& A) { foo_impl(A); } |
|||
* \endcode |
|||
* |
|||
* |
|||
* \sa PlainObjectBase::Map(), \ref TopicStorageOrders |
|||
*/ |
|||
|
|||
namespace internal { |
|||
|
|||
template<typename _PlainObjectType, int _Options, typename _StrideType> |
|||
struct traits<Ref<_PlainObjectType, _Options, _StrideType> > |
|||
: public traits<Map<_PlainObjectType, _Options, _StrideType> > |
|||
{ |
|||
typedef _PlainObjectType PlainObjectType; |
|||
typedef _StrideType StrideType; |
|||
enum { |
|||
Options = _Options |
|||
}; |
|||
|
|||
template<typename Derived> struct match { |
|||
enum { |
|||
HasDirectAccess = internal::has_direct_access<Derived>::ret, |
|||
StorageOrderMatch = PlainObjectType::IsVectorAtCompileTime || ((PlainObjectType::Flags&RowMajorBit)==(Derived::Flags&RowMajorBit)), |
|||
InnerStrideMatch = int(StrideType::InnerStrideAtCompileTime)==int(Dynamic) |
|||
|| int(StrideType::InnerStrideAtCompileTime)==int(Derived::InnerStrideAtCompileTime) |
|||
|| (int(StrideType::InnerStrideAtCompileTime)==0 && int(Derived::InnerStrideAtCompileTime)==1), |
|||
OuterStrideMatch = Derived::IsVectorAtCompileTime |
|||
|| int(StrideType::OuterStrideAtCompileTime)==int(Dynamic) || int(StrideType::OuterStrideAtCompileTime)==int(Derived::OuterStrideAtCompileTime), |
|||
AlignmentMatch = (_Options!=Aligned) || ((PlainObjectType::Flags&AlignedBit)==0) || ((traits<Derived>::Flags&AlignedBit)==AlignedBit), |
|||
MatchAtCompileTime = HasDirectAccess && StorageOrderMatch && InnerStrideMatch && OuterStrideMatch && AlignmentMatch |
|||
}; |
|||
typedef typename internal::conditional<MatchAtCompileTime,internal::true_type,internal::false_type>::type type; |
|||
}; |
|||
|
|||
}; |
|||
|
|||
template<typename Derived> |
|||
struct traits<RefBase<Derived> > : public traits<Derived> {}; |
|||
|
|||
} |
|||
|
|||
template<typename Derived> class RefBase |
|||
: public MapBase<Derived> |
|||
{ |
|||
typedef typename internal::traits<Derived>::PlainObjectType PlainObjectType; |
|||
typedef typename internal::traits<Derived>::StrideType StrideType; |
|||
|
|||
public: |
|||
|
|||
typedef MapBase<Derived> Base; |
|||
EIGEN_DENSE_PUBLIC_INTERFACE(RefBase) |
|||
|
|||
inline Index innerStride() const |
|||
{ |
|||
return StrideType::InnerStrideAtCompileTime != 0 ? m_stride.inner() : 1; |
|||
} |
|||
|
|||
inline Index outerStride() const |
|||
{ |
|||
return StrideType::OuterStrideAtCompileTime != 0 ? m_stride.outer() |
|||
: IsVectorAtCompileTime ? this->size() |
|||
: int(Flags)&RowMajorBit ? this->cols() |
|||
: this->rows(); |
|||
} |
|||
|
|||
RefBase() |
|||
: Base(0,RowsAtCompileTime==Dynamic?0:RowsAtCompileTime,ColsAtCompileTime==Dynamic?0:ColsAtCompileTime), |
|||
// Stride<> does not allow default ctor for Dynamic strides, so let' initialize it with dummy values: |
|||
m_stride(StrideType::OuterStrideAtCompileTime==Dynamic?0:StrideType::OuterStrideAtCompileTime, |
|||
StrideType::InnerStrideAtCompileTime==Dynamic?0:StrideType::InnerStrideAtCompileTime) |
|||
{} |
|||
|
|||
protected: |
|||
|
|||
typedef Stride<StrideType::OuterStrideAtCompileTime,StrideType::InnerStrideAtCompileTime> StrideBase; |
|||
|
|||
template<typename Expression> |
|||
void construct(Expression& expr) |
|||
{ |
|||
if(PlainObjectType::RowsAtCompileTime==1) |
|||
{ |
|||
eigen_assert(expr.rows()==1 || expr.cols()==1); |
|||
::new (static_cast<Base*>(this)) Base(expr.data(), 1, expr.size()); |
|||
} |
|||
else if(PlainObjectType::ColsAtCompileTime==1) |
|||
{ |
|||
eigen_assert(expr.rows()==1 || expr.cols()==1); |
|||
::new (static_cast<Base*>(this)) Base(expr.data(), expr.size(), 1); |
|||
} |
|||
else |
|||
::new (static_cast<Base*>(this)) Base(expr.data(), expr.rows(), expr.cols()); |
|||
::new (&m_stride) StrideBase(StrideType::OuterStrideAtCompileTime==0?0:expr.outerStride(), |
|||
StrideType::InnerStrideAtCompileTime==0?0:expr.innerStride()); |
|||
} |
|||
|
|||
StrideBase m_stride; |
|||
}; |
|||
|
|||
|
|||
template<typename PlainObjectType, int Options, typename StrideType> class Ref |
|||
: public RefBase<Ref<PlainObjectType, Options, StrideType> > |
|||
{ |
|||
typedef internal::traits<Ref> Traits; |
|||
public: |
|||
|
|||
typedef RefBase<Ref> Base; |
|||
EIGEN_DENSE_PUBLIC_INTERFACE(Ref) |
|||
|
|||
|
|||
#ifndef EIGEN_PARSED_BY_DOXYGEN |
|||
template<typename Derived> |
|||
inline Ref(PlainObjectBase<Derived>& expr, |
|||
typename internal::enable_if<bool(Traits::template match<Derived>::MatchAtCompileTime),Derived>::type* = 0) |
|||
{ |
|||
Base::construct(expr); |
|||
} |
|||
template<typename Derived> |
|||
inline Ref(const DenseBase<Derived>& expr, |
|||
typename internal::enable_if<bool(internal::is_lvalue<Derived>::value&&bool(Traits::template match<Derived>::MatchAtCompileTime)),Derived>::type* = 0, |
|||
int = Derived::ThisConstantIsPrivateInPlainObjectBase) |
|||
#else |
|||
template<typename Derived> |
|||
inline Ref(DenseBase<Derived>& expr) |
|||
#endif |
|||
{ |
|||
Base::construct(expr.const_cast_derived()); |
|||
} |
|||
|
|||
EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Ref) |
|||
|
|||
}; |
|||
|
|||
// this is the const ref version |
|||
template<typename PlainObjectType, int Options, typename StrideType> class Ref<const PlainObjectType, Options, StrideType> |
|||
: public RefBase<Ref<const PlainObjectType, Options, StrideType> > |
|||
{ |
|||
typedef internal::traits<Ref> Traits; |
|||
public: |
|||
|
|||
typedef RefBase<Ref> Base; |
|||
EIGEN_DENSE_PUBLIC_INTERFACE(Ref) |
|||
|
|||
template<typename Derived> |
|||
inline Ref(const DenseBase<Derived>& expr) |
|||
{ |
|||
// std::cout << match_helper<Derived>::HasDirectAccess << "," << match_helper<Derived>::OuterStrideMatch << "," << match_helper<Derived>::InnerStrideMatch << "\n"; |
|||
// std::cout << int(StrideType::OuterStrideAtCompileTime) << " - " << int(Derived::OuterStrideAtCompileTime) << "\n"; |
|||
// std::cout << int(StrideType::InnerStrideAtCompileTime) << " - " << int(Derived::InnerStrideAtCompileTime) << "\n"; |
|||
construct(expr.derived(), typename Traits::template match<Derived>::type()); |
|||
} |
|||
|
|||
protected: |
|||
|
|||
template<typename Expression> |
|||
void construct(const Expression& expr,internal::true_type) |
|||
{ |
|||
Base::construct(expr); |
|||
} |
|||
|
|||
template<typename Expression> |
|||
void construct(const Expression& expr, internal::false_type) |
|||
{ |
|||
// std::cout << "Ref: copy\n"; |
|||
m_object = expr; |
|||
Base::construct(m_object); |
|||
} |
|||
|
|||
protected: |
|||
PlainObjectType m_object; |
|||
}; |
|||
|
|||
} // end namespace Eigen |
|||
|
|||
#endif // EIGEN_REF_H |
@ -0,0 +1,339 @@ |
|||
// This file is part of Eigen, a lightweight C++ template library |
|||
// for linear algebra. |
|||
// |
|||
// Copyright (C) 2012 Gael Guennebaud <gael.guennebaud@inria.fr> |
|||
// Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk> |
|||
// |
|||
// This Source Code Form is subject to the terms of the Mozilla |
|||
// Public License v. 2.0. If a copy of the MPL was not distributed |
|||
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. |
|||
|
|||
#ifndef EIGEN_GENERALIZEDEIGENSOLVER_H |
|||
#define EIGEN_GENERALIZEDEIGENSOLVER_H |
|||
|
|||
#include "./RealQZ.h" |
|||
|
|||
namespace Eigen { |
|||
|
|||
/** \eigenvalues_module \ingroup Eigenvalues_Module |
|||
* |
|||
* |
|||
* \class GeneralizedEigenSolver |
|||
* |
|||
* \brief Computes the generalized eigenvalues and eigenvectors of a pair of general matrices |
|||
* |
|||
* \tparam _MatrixType the type of the matrices of which we are computing the |
|||
* eigen-decomposition; this is expected to be an instantiation of the Matrix |
|||
* class template. Currently, only real matrices are supported. |
|||
* |
|||
* The generalized eigenvalues and eigenvectors of a matrix pair \f$ A \f$ and \f$ B \f$ are scalars |
|||
* \f$ \lambda \f$ and vectors \f$ v \f$ such that \f$ Av = \lambda Bv \f$. If |
|||
* \f$ D \f$ is a diagonal matrix with the eigenvalues on the diagonal, and |
|||
* \f$ V \f$ is a matrix with the eigenvectors as its columns, then \f$ A V = |
|||
* B V D \f$. The matrix \f$ V \f$ is almost always invertible, in which case we |
|||
* have \f$ A = B V D V^{-1} \f$. This is called the generalized eigen-decomposition. |
|||
* |
|||
* The generalized eigenvalues and eigenvectors of a matrix pair may be complex, even when the |
|||
* matrices are real. Moreover, the generalized eigenvalue might be infinite if the matrix B is |
|||
* singular. To workaround this difficulty, the eigenvalues are provided as a pair of complex \f$ \alpha \f$ |
|||
* and real \f$ \beta \f$ such that: \f$ \lambda_i = \alpha_i / \beta_i \f$. If \f$ \beta_i \f$ is (nearly) zero, |
|||
* then one can consider the well defined left eigenvalue \f$ \mu = \beta_i / \alpha_i\f$ such that: |
|||
* \f$ \mu_i A v_i = B v_i \f$, or even \f$ \mu_i u_i^T A = u_i^T B \f$ where \f$ u_i \f$ is |
|||
* called the left eigenvector. |
|||
* |
|||
* Call the function compute() to compute the generalized eigenvalues and eigenvectors of |
|||
* a given matrix pair. Alternatively, you can use the |
|||
* GeneralizedEigenSolver(const MatrixType&, const MatrixType&, bool) constructor which computes the |
|||
* eigenvalues and eigenvectors at construction time. Once the eigenvalue and |
|||
* eigenvectors are computed, they can be retrieved with the eigenvalues() and |
|||
* eigenvectors() functions. |
|||
* |
|||
* Here is an usage example of this class: |
|||
* Example: \include GeneralizedEigenSolver.cpp |
|||
* Output: \verbinclude GeneralizedEigenSolver.out |
|||
* |
|||
* \sa MatrixBase::eigenvalues(), class ComplexEigenSolver, class SelfAdjointEigenSolver |
|||
*/ |
|||
template<typename _MatrixType> class GeneralizedEigenSolver |
|||
{ |
|||
public: |
|||
|
|||
/** \brief Synonym for the template parameter \p _MatrixType. */ |
|||
typedef _MatrixType MatrixType; |
|||
|
|||
enum { |
|||
RowsAtCompileTime = MatrixType::RowsAtCompileTime, |
|||
ColsAtCompileTime = MatrixType::ColsAtCompileTime, |
|||
Options = MatrixType::Options, |
|||
MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, |
|||
MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime |
|||
}; |
|||
|
|||
/** \brief Scalar type for matrices of type #MatrixType. */ |
|||
typedef typename MatrixType::Scalar Scalar; |
|||
typedef typename NumTraits<Scalar>::Real RealScalar; |
|||
typedef typename MatrixType::Index Index; |
|||
|
|||
/** \brief Complex scalar type for #MatrixType. |
|||
* |
|||
* This is \c std::complex<Scalar> if #Scalar is real (e.g., |
|||
* \c float or \c double) and just \c Scalar if #Scalar is |
|||
* complex. |
|||
*/ |
|||
typedef std::complex<RealScalar> ComplexScalar; |
|||
|
|||
/** \brief Type for vector of real scalar values eigenvalues as returned by betas(). |
|||
* |
|||
* This is a column vector with entries of type #Scalar. |
|||
* The length of the vector is the size of #MatrixType. |
|||
*/ |
|||
typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> VectorType; |
|||
|
|||
/** \brief Type for vector of complex scalar values eigenvalues as returned by betas(). |
|||
* |
|||
* This is a column vector with entries of type #ComplexScalar. |
|||
* The length of the vector is the size of #MatrixType. |
|||
*/ |
|||
typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ComplexVectorType; |
|||
|
|||
/** \brief Expression type for the eigenvalues as returned by eigenvalues(). |
|||
*/ |
|||
typedef CwiseBinaryOp<internal::scalar_quotient_op<ComplexScalar,Scalar>,ComplexVectorType,VectorType> EigenvalueType; |
|||
|
|||
/** \brief Type for matrix of eigenvectors as returned by eigenvectors(). |
|||
* |
|||
* This is a square matrix with entries of type #ComplexScalar. |
|||
* The size is the same as the size of #MatrixType. |
|||
*/ |
|||
typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorsType; |
|||
|
|||
/** \brief Default constructor. |
|||
* |
|||
* The default constructor is useful in cases in which the user intends to |
|||
* perform decompositions via EigenSolver::compute(const MatrixType&, bool). |
|||
* |
|||
* \sa compute() for an example. |
|||
*/ |
|||
GeneralizedEigenSolver() : m_eivec(), m_alphas(), m_betas(), m_isInitialized(false), m_realQZ(), m_matS(), m_tmp() {} |
|||
|
|||
/** \brief Default constructor with memory preallocation |
|||
* |
|||
* Like the default constructor but with preallocation of the internal data |
|||
* according to the specified problem \a size. |
|||
* \sa GeneralizedEigenSolver() |
|||
*/ |
|||
GeneralizedEigenSolver(Index size) |
|||
: m_eivec(size, size), |
|||
m_alphas(size), |
|||
m_betas(size), |
|||
m_isInitialized(false), |
|||
m_eigenvectorsOk(false), |
|||
m_realQZ(size), |
|||
m_matS(size, size), |
|||
m_tmp(size) |
|||
{} |
|||
|
|||
/** \brief Constructor; computes the generalized eigendecomposition of given matrix pair. |
|||
* |
|||
* \param[in] A Square matrix whose eigendecomposition is to be computed. |
|||
* \param[in] B Square matrix whose eigendecomposition is to be computed. |
|||
* \param[in] computeEigenvectors If true, both the eigenvectors and the |
|||
* eigenvalues are computed; if false, only the eigenvalues are computed. |
|||
* |
|||
* This constructor calls compute() to compute the generalized eigenvalues |
|||
* and eigenvectors. |
|||
* |
|||
* \sa compute() |
|||
*/ |
|||
GeneralizedEigenSolver(const MatrixType& A, const MatrixType& B, bool computeEigenvectors = true) |
|||
: m_eivec(A.rows(), A.cols()), |
|||
m_alphas(A.cols()), |
|||
m_betas(A.cols()), |
|||
m_isInitialized(false), |
|||
m_eigenvectorsOk(false), |
|||
m_realQZ(A.cols()), |
|||
m_matS(A.rows(), A.cols()), |
|||
m_tmp(A.cols()) |
|||
{ |
|||
compute(A, B, computeEigenvectors); |
|||
} |
|||
|
|||
/* \brief Returns the computed generalized eigenvectors. |
|||
* |
|||
* \returns %Matrix whose columns are the (possibly complex) eigenvectors. |
|||
* |
|||
* \pre Either the constructor |
|||
* GeneralizedEigenSolver(const MatrixType&,const MatrixType&, bool) or the member function |
|||
* compute(const MatrixType&, const MatrixType& bool) has been called before, and |
|||
* \p computeEigenvectors was set to true (the default). |
|||
* |
|||
* Column \f$ k \f$ of the returned matrix is an eigenvector corresponding |
|||
* to eigenvalue number \f$ k \f$ as returned by eigenvalues(). The |
|||
* eigenvectors are normalized to have (Euclidean) norm equal to one. The |
|||
* matrix returned by this function is the matrix \f$ V \f$ in the |
|||
* generalized eigendecomposition \f$ A = B V D V^{-1} \f$, if it exists. |
|||
* |
|||
* \sa eigenvalues() |
|||
*/ |
|||
// EigenvectorsType eigenvectors() const; |
|||
|
|||
/** \brief Returns an expression of the computed generalized eigenvalues. |
|||
* |
|||
* \returns An expression of the column vector containing the eigenvalues. |
|||
* |
|||
* It is a shortcut for \code this->alphas().cwiseQuotient(this->betas()); \endcode |
|||
* Not that betas might contain zeros. It is therefore not recommended to use this function, |
|||
* but rather directly deal with the alphas and betas vectors. |
|||
* |
|||
* \pre Either the constructor |
|||
* GeneralizedEigenSolver(const MatrixType&,const MatrixType&,bool) or the member function |
|||
* compute(const MatrixType&,const MatrixType&,bool) has been called before. |
|||
* |
|||
* The eigenvalues are repeated according to their algebraic multiplicity, |
|||
* so there are as many eigenvalues as rows in the matrix. The eigenvalues |
|||
* are not sorted in any particular order. |
|||
* |
|||
* \sa alphas(), betas(), eigenvectors() |
|||
*/ |
|||
EigenvalueType eigenvalues() const |
|||
{ |
|||
eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized."); |
|||
return EigenvalueType(m_alphas,m_betas); |
|||
} |
|||
|
|||
/** \returns A const reference to the vectors containing the alpha values |
|||
* |
|||
* This vector permits to reconstruct the j-th eigenvalues as alphas(i)/betas(j). |
|||
* |
|||
* \sa betas(), eigenvalues() */ |
|||
ComplexVectorType alphas() const |
|||
{ |
|||
eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized."); |
|||
return m_alphas; |
|||
} |
|||
|
|||
/** \returns A const reference to the vectors containing the beta values |
|||
* |
|||
* This vector permits to reconstruct the j-th eigenvalues as alphas(i)/betas(j). |
|||
* |
|||
* \sa alphas(), eigenvalues() */ |
|||
VectorType betas() const |
|||
{ |
|||
eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized."); |
|||
return m_betas; |
|||
} |
|||
|
|||
/** \brief Computes generalized eigendecomposition of given matrix. |
|||
* |
|||
* \param[in] A Square matrix whose eigendecomposition is to be computed. |
|||
* \param[in] B Square matrix whose eigendecomposition is to be computed. |
|||
* \param[in] computeEigenvectors If true, both the eigenvectors and the |
|||
* eigenvalues are computed; if false, only the eigenvalues are |
|||
* computed. |
|||
* \returns Reference to \c *this |
|||
* |
|||
* This function computes the eigenvalues of the real matrix \p matrix. |
|||
* The eigenvalues() function can be used to retrieve them. If |
|||
* \p computeEigenvectors is true, then the eigenvectors are also computed |
|||
* and can be retrieved by calling eigenvectors(). |
|||
* |
|||
* The matrix is first reduced to real generalized Schur form using the RealQZ |
|||
* class. The generalized Schur decomposition is then used to compute the eigenvalues |
|||
* and eigenvectors. |
|||
* |
|||
* The cost of the computation is dominated by the cost of the |
|||
* generalized Schur decomposition. |
|||
* |
|||
* This method reuses of the allocated data in the GeneralizedEigenSolver object. |
|||
*/ |
|||
GeneralizedEigenSolver& compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors = true); |
|||
|
|||
ComputationInfo info() const |
|||
{ |
|||
eigen_assert(m_isInitialized && "EigenSolver is not initialized."); |
|||
return m_realQZ.info(); |
|||
} |
|||
|
|||
/** Sets the maximal number of iterations allowed. |
|||
*/ |
|||
GeneralizedEigenSolver& setMaxIterations(Index maxIters) |
|||
{ |
|||
m_realQZ.setMaxIterations(maxIters); |
|||
return *this; |
|||
} |
|||
|
|||
protected: |
|||
MatrixType m_eivec; |
|||
ComplexVectorType m_alphas; |
|||
VectorType m_betas; |
|||
bool m_isInitialized; |
|||
bool m_eigenvectorsOk; |
|||
RealQZ<MatrixType> m_realQZ; |
|||
MatrixType m_matS; |
|||
|
|||
typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType; |
|||
ColumnVectorType m_tmp; |
|||
}; |
|||
|
|||
//template<typename MatrixType> |
|||
//typename GeneralizedEigenSolver<MatrixType>::EigenvectorsType GeneralizedEigenSolver<MatrixType>::eigenvectors() const |
|||
//{ |
|||
// eigen_assert(m_isInitialized && "EigenSolver is not initialized."); |
|||
// eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); |
|||
// Index n = m_eivec.cols(); |
|||
// EigenvectorsType matV(n,n); |
|||
// // TODO |
|||
// return matV; |
|||
//} |
|||
|
|||
template<typename MatrixType> |
|||
GeneralizedEigenSolver<MatrixType>& |
|||
GeneralizedEigenSolver<MatrixType>::compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors) |
|||
{ |
|||
eigen_assert(A.cols() == A.rows() && B.cols() == A.rows() && B.cols() == B.rows()); |
|||
|
|||
// Reduce to generalized real Schur form: |
|||
// A = Q S Z and B = Q T Z |
|||
m_realQZ.compute(A, B, computeEigenvectors); |
|||
|
|||
if (m_realQZ.info() == Success) |
|||
{ |
|||
m_matS = m_realQZ.matrixS(); |
|||
if (computeEigenvectors) |
|||
m_eivec = m_realQZ.matrixZ().transpose(); |
|||
|
|||
// Compute eigenvalues from matS |
|||
m_alphas.resize(A.cols()); |
|||
m_betas.resize(A.cols()); |
|||
Index i = 0; |
|||
while (i < A.cols()) |
|||
{ |
|||
if (i == A.cols() - 1 || m_matS.coeff(i+1, i) == Scalar(0)) |
|||
{ |
|||
m_alphas.coeffRef(i) = m_matS.coeff(i, i); |
|||
m_betas.coeffRef(i) = m_realQZ.matrixT().coeff(i,i); |
|||
++i; |
|||
} |
|||
else |
|||
{ |
|||
Scalar p = Scalar(0.5) * (m_matS.coeff(i, i) - m_matS.coeff(i+1, i+1)); |
|||
Scalar z = internal::sqrt(internal::abs(p * p + m_matS.coeff(i+1, i) * m_matS.coeff(i, i+1))); |
|||
m_alphas.coeffRef(i) = ComplexScalar(m_matS.coeff(i+1, i+1) + p, z); |
|||
m_alphas.coeffRef(i+1) = ComplexScalar(m_matS.coeff(i+1, i+1) + p, -z); |
|||
|
|||
m_betas.coeffRef(i) = m_realQZ.matrixT().coeff(i,i); |
|||
m_betas.coeffRef(i+1) = m_realQZ.matrixT().coeff(i,i); |
|||
i += 2; |
|||
} |
|||
} |
|||
} |
|||
|
|||
m_isInitialized = true; |
|||
m_eigenvectorsOk = false;//computeEigenvectors; |
|||
|
|||
return *this; |
|||
} |
|||
|
|||
} // end namespace Eigen |
|||
|
|||
#endif // EIGEN_GENERALIZEDEIGENSOLVER_H |
@ -0,0 +1,618 @@ |
|||
// This file is part of Eigen, a lightweight C++ template library |
|||
// for linear algebra. |
|||
// |
|||
// Copyright (C) 2012 Alexey Korepanov <kaikaikai@yandex.ru> |
|||
// |
|||
// This Source Code Form is subject to the terms of the Mozilla |
|||
// Public License v. 2.0. If a copy of the MPL was not distributed |
|||
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. |
|||
|
|||
#ifndef EIGEN_REAL_QZ_H |
|||
#define EIGEN_REAL_QZ_H |
|||
|
|||
namespace Eigen { |
|||
|
|||
/** \eigenvalues_module \ingroup Eigenvalues_Module |
|||
* |
|||
* |
|||
* \class RealQZ |
|||
* |
|||
* \brief Performs a real QZ decomposition of a pair of square matrices |
|||
* |
|||
* \tparam _MatrixType the type of the matrix of which we are computing the |
|||
* real QZ decomposition; this is expected to be an instantiation of the |
|||
* Matrix class template. |
|||
* |
|||
* Given a real square matrices A and B, this class computes the real QZ |
|||
* decomposition: \f$ A = Q S Z \f$, \f$ B = Q T Z \f$ where Q and Z are |
|||
* real orthogonal matrixes, T is upper-triangular matrix, and S is upper |
|||
* quasi-triangular matrix. An orthogonal matrix is a matrix whose |
|||
* inverse is equal to its transpose, \f$ U^{-1} = U^T \f$. A quasi-triangular |
|||
* matrix is a block-triangular matrix whose diagonal consists of 1-by-1 |
|||
* blocks and 2-by-2 blocks where further reduction is impossible due to |
|||
* complex eigenvalues. |
|||
* |
|||
* The eigenvalues of the pencil \f$ A - z B \f$ can be obtained from |
|||
* 1x1 and 2x2 blocks on the diagonals of S and T. |
|||
* |
|||
* Call the function compute() to compute the real QZ decomposition of a |
|||
* given pair of matrices. Alternatively, you can use the |
|||
* RealQZ(const MatrixType& B, const MatrixType& B, bool computeQZ) |
|||
* constructor which computes the real QZ decomposition at construction |
|||
* time. Once the decomposition is computed, you can use the matrixS(), |
|||
* matrixT(), matrixQ() and matrixZ() functions to retrieve the matrices |
|||
* S, T, Q and Z in the decomposition. If computeQZ==false, some time |
|||
* is saved by not computing matrices Q and Z. |
|||
* |
|||
* Example: \include RealQZ_compute.cpp |
|||
* Output: \include RealQZ_compute.out |
|||
* |
|||
* \note The implementation is based on the algorithm in "Matrix Computations" |
|||
* by Gene H. Golub and Charles F. Van Loan, and a paper "An algorithm for |
|||
* generalized eigenvalue problems" by C.B.Moler and G.W.Stewart. |
|||
* |
|||
* \sa class RealSchur, class ComplexSchur, class EigenSolver, class ComplexEigenSolver |
|||
*/ |
|||
|
|||
template<typename _MatrixType> class RealQZ |
|||
{ |
|||
public: |
|||
typedef _MatrixType MatrixType; |
|||
enum { |
|||
RowsAtCompileTime = MatrixType::RowsAtCompileTime, |
|||
ColsAtCompileTime = MatrixType::ColsAtCompileTime, |
|||
Options = MatrixType::Options, |
|||
MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, |
|||
MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime |
|||
}; |
|||
typedef typename MatrixType::Scalar Scalar; |
|||
typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar; |
|||
typedef typename MatrixType::Index Index; |
|||
|
|||
typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> EigenvalueType; |
|||
typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType; |
|||
|
|||
/** \brief Default constructor. |
|||
* |
|||
* \param [in] size Positive integer, size of the matrix whose QZ decomposition will be computed. |
|||
* |
|||
* The default constructor is useful in cases in which the user intends to |
|||
* perform decompositions via compute(). The \p size parameter is only |
|||
* used as a hint. It is not an error to give a wrong \p size, but it may |
|||
* impair performance. |
|||
* |
|||
* \sa compute() for an example. |
|||
*/ |
|||
RealQZ(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime) : |
|||
m_S(size, size), |
|||
m_T(size, size), |
|||
m_Q(size, size), |
|||
m_Z(size, size), |
|||
m_workspace(size*2), |
|||
m_maxIters(400), |
|||
m_isInitialized(false) |
|||
{ } |
|||
|
|||
/** \brief Constructor; computes real QZ decomposition of given matrices |
|||
* |
|||
* \param[in] A Matrix A. |
|||
* \param[in] B Matrix B. |
|||
* \param[in] computeQZ If false, A and Z are not computed. |
|||
* |
|||
* This constructor calls compute() to compute the QZ decomposition. |
|||
*/ |
|||
RealQZ(const MatrixType& A, const MatrixType& B, bool computeQZ = true) : |
|||
m_S(A.rows(),A.cols()), |
|||
m_T(A.rows(),A.cols()), |
|||
m_Q(A.rows(),A.cols()), |
|||
m_Z(A.rows(),A.cols()), |
|||
m_workspace(A.rows()*2), |
|||
m_maxIters(400), |
|||
m_isInitialized(false) { |
|||
compute(A, B, computeQZ); |
|||
} |
|||
|
|||
/** \brief Returns matrix Q in the QZ decomposition. |
|||
* |
|||
* \returns A const reference to the matrix Q. |
|||
*/ |
|||
const MatrixType& matrixQ() const { |
|||
eigen_assert(m_isInitialized && "RealQZ is not initialized."); |
|||
eigen_assert(m_computeQZ && "The matrices Q and Z have not been computed during the QZ decomposition."); |
|||
return m_Q; |
|||
} |
|||
|
|||
/** \brief Returns matrix Z in the QZ decomposition. |
|||
* |
|||
* \returns A const reference to the matrix Z. |
|||
*/ |
|||
const MatrixType& matrixZ() const { |
|||
eigen_assert(m_isInitialized && "RealQZ is not initialized."); |
|||
eigen_assert(m_computeQZ && "The matrices Q and Z have not been computed during the QZ decomposition."); |
|||
return m_Z; |
|||
} |
|||
|
|||
/** \brief Returns matrix S in the QZ decomposition. |
|||
* |
|||
* \returns A const reference to the matrix S. |
|||
*/ |
|||
const MatrixType& matrixS() const { |
|||
eigen_assert(m_isInitialized && "RealQZ is not initialized."); |
|||
return m_S; |
|||
} |
|||
|
|||
/** \brief Returns matrix S in the QZ decomposition. |
|||
* |
|||
* \returns A const reference to the matrix S. |
|||
*/ |
|||
const MatrixType& matrixT() const { |
|||
eigen_assert(m_isInitialized && "RealQZ is not initialized."); |
|||
return m_T; |
|||
} |
|||
|
|||
/** \brief Computes QZ decomposition of given matrix. |
|||
* |
|||
* \param[in] A Matrix A. |
|||
* \param[in] B Matrix B. |
|||
* \param[in] computeQZ If false, A and Z are not computed. |
|||
* \returns Reference to \c *this |
|||
*/ |
|||
RealQZ& compute(const MatrixType& A, const MatrixType& B, bool computeQZ = true); |
|||
|
|||
/** \brief Reports whether previous computation was successful. |
|||
* |
|||
* \returns \c Success if computation was succesful, \c NoConvergence otherwise. |
|||
*/ |
|||
ComputationInfo info() const |
|||
{ |
|||
eigen_assert(m_isInitialized && "RealQZ is not initialized."); |
|||
return m_info; |
|||
} |
|||
|
|||
/** \brief Returns number of performed QR-like iterations. |
|||
*/ |
|||
Index iterations() const |
|||
{ |
|||
eigen_assert(m_isInitialized && "RealQZ is not initialized."); |
|||
return m_global_iter; |
|||
} |
|||
|
|||
/** Sets the maximal number of iterations allowed to converge to one eigenvalue |
|||
* or decouple the problem. |
|||
*/ |
|||
RealQZ& setMaxIterations(Index maxIters) |
|||
{ |
|||
m_maxIters = maxIters; |
|||
return *this; |
|||
} |
|||
|
|||
private: |
|||
|
|||
MatrixType m_S, m_T, m_Q, m_Z; |
|||
Matrix<Scalar,Dynamic,1> m_workspace; |
|||
ComputationInfo m_info; |
|||
Index m_maxIters; |
|||
bool m_isInitialized; |
|||
bool m_computeQZ; |
|||
Scalar m_normOfT, m_normOfS; |
|||
Index m_global_iter; |
|||
|
|||
typedef Matrix<Scalar,3,1> Vector3s; |
|||
typedef Matrix<Scalar,2,1> Vector2s; |
|||
typedef Matrix<Scalar,2,2> Matrix2s; |
|||
typedef JacobiRotation<Scalar> JRs; |
|||
|
|||
void hessenbergTriangular(); |
|||
void computeNorms(); |
|||
Index findSmallSubdiagEntry(Index iu); |
|||
Index findSmallDiagEntry(Index f, Index l); |
|||
void splitOffTwoRows(Index i); |
|||
void pushDownZero(Index z, Index f, Index l); |
|||
void step(Index f, Index l, Index iter); |
|||
|
|||
}; // RealQZ |
|||
|
|||
/** \internal Reduces S and T to upper Hessenberg - triangular form */ |
|||
template<typename MatrixType> |
|||
void RealQZ<MatrixType>::hessenbergTriangular() |
|||
{ |
|||
|
|||
const Index dim = m_S.cols(); |
|||
|
|||
// perform QR decomposition of T, overwrite T with R, save Q |
|||
HouseholderQR<MatrixType> qrT(m_T); |
|||
m_T = qrT.matrixQR(); |
|||
m_T.template triangularView<StrictlyLower>().setZero(); |
|||
m_Q = qrT.householderQ(); |
|||
// overwrite S with Q* S |
|||
m_S.applyOnTheLeft(m_Q.adjoint()); |
|||
// init Z as Identity |
|||
if (m_computeQZ) |
|||
m_Z = MatrixType::Identity(dim,dim); |
|||
// reduce S to upper Hessenberg with Givens rotations |
|||
for (Index j=0; j<=dim-3; j++) { |
|||
for (Index i=dim-1; i>=j+2; i--) { |
|||
JRs G; |
|||
// kill S(i,j) |
|||
if(m_S.coeff(i,j) != 0) |
|||
{ |
|||
G.makeGivens(m_S.coeff(i-1,j), m_S.coeff(i,j), &m_S.coeffRef(i-1, j)); |
|||
m_S.coeffRef(i,j) = Scalar(0.0); |
|||
m_S.rightCols(dim-j-1).applyOnTheLeft(i-1,i,G.adjoint()); |
|||
m_T.rightCols(dim-i+1).applyOnTheLeft(i-1,i,G.adjoint()); |
|||
} |
|||
// update Q |
|||
if (m_computeQZ) |
|||
m_Q.applyOnTheRight(i-1,i,G); |
|||
// kill T(i,i-1) |
|||
if(m_T.coeff(i,i-1)!=Scalar(0)) |
|||
{ |
|||
G.makeGivens(m_T.coeff(i,i), m_T.coeff(i,i-1), &m_T.coeffRef(i,i)); |
|||
m_T.coeffRef(i,i-1) = Scalar(0.0); |
|||
m_S.applyOnTheRight(i,i-1,G); |
|||
m_T.topRows(i).applyOnTheRight(i,i-1,G); |
|||
} |
|||
// update Z |
|||
if (m_computeQZ) |
|||
m_Z.applyOnTheLeft(i,i-1,G.adjoint()); |
|||
} |
|||
} |
|||
} |
|||
|
|||
/** \internal Computes vector L1 norms of S and T when in Hessenberg-Triangular form already */ |
|||
template<typename MatrixType> |
|||
inline void RealQZ<MatrixType>::computeNorms() |
|||
{ |
|||
const Index size = m_S.cols(); |
|||
m_normOfS = Scalar(0.0); |
|||
m_normOfT = Scalar(0.0); |
|||
for (Index j = 0; j < size; ++j) |
|||
{ |
|||
m_normOfS += m_S.col(j).segment(0, (std::min)(size,j+2)).cwiseAbs().sum(); |
|||
m_normOfT += m_T.row(j).segment(j, size - j).cwiseAbs().sum(); |
|||
} |
|||
} |
|||
|
|||
|
|||
/** \internal Look for single small sub-diagonal element S(res, res-1) and return res (or 0) */ |
|||
template<typename MatrixType> |
|||
inline typename MatrixType::Index RealQZ<MatrixType>::findSmallSubdiagEntry(Index iu) |
|||
{ |
|||
Index res = iu; |
|||
while (res > 0) |
|||
{ |
|||
Scalar s = internal::abs(m_S.coeff(res-1,res-1)) + internal::abs(m_S.coeff(res,res)); |
|||
if (s == Scalar(0.0)) |
|||
s = m_normOfS; |
|||
if (internal::abs(m_S.coeff(res,res-1)) < NumTraits<Scalar>::epsilon() * s) |
|||
break; |
|||
res--; |
|||
} |
|||
return res; |
|||
} |
|||
|
|||
/** \internal Look for single small diagonal element T(res, res) for res between f and l, and return res (or f-1) */ |
|||
template<typename MatrixType> |
|||
inline typename MatrixType::Index RealQZ<MatrixType>::findSmallDiagEntry(Index f, Index l) |
|||
{ |
|||
Index res = l; |
|||
while (res >= f) { |
|||
if (internal::abs(m_T.coeff(res,res)) <= NumTraits<Scalar>::epsilon() * m_normOfT) |
|||
break; |
|||
res--; |
|||
} |
|||
return res; |
|||
} |
|||
|
|||
/** \internal decouple 2x2 diagonal block in rows i, i+1 if eigenvalues are real */ |
|||
template<typename MatrixType> |
|||
inline void RealQZ<MatrixType>::splitOffTwoRows(Index i) |
|||
{ |
|||
const Index dim=m_S.cols(); |
|||
if (internal::abs(m_S.coeff(i+1,i)==Scalar(0))) |
|||
return; |
|||
Index z = findSmallDiagEntry(i,i+1); |
|||
if (z==i-1) |
|||
{ |
|||
// block of (S T^{-1}) |
|||
Matrix2s STi = m_T.template block<2,2>(i,i).template triangularView<Upper>(). |
|||
template solve<OnTheRight>(m_S.template block<2,2>(i,i)); |
|||
Scalar p = Scalar(0.5)*(STi(0,0)-STi(1,1)); |
|||
Scalar q = p*p + STi(1,0)*STi(0,1); |
|||
if (q>=0) { |
|||
Scalar z = internal::sqrt(q); |
|||
// one QR-like iteration for ABi - lambda I |
|||
// is enough - when we know exact eigenvalue in advance, |
|||
// convergence is immediate |
|||
JRs G; |
|||
if (p>=0) |
|||
G.makeGivens(p + z, STi(1,0)); |
|||
else |
|||
G.makeGivens(p - z, STi(1,0)); |
|||
m_S.rightCols(dim-i).applyOnTheLeft(i,i+1,G.adjoint()); |
|||
m_T.rightCols(dim-i).applyOnTheLeft(i,i+1,G.adjoint()); |
|||
// update Q |
|||
if (m_computeQZ) |
|||
m_Q.applyOnTheRight(i,i+1,G); |
|||
|
|||
G.makeGivens(m_T.coeff(i+1,i+1), m_T.coeff(i+1,i)); |
|||
m_S.topRows(i+2).applyOnTheRight(i+1,i,G); |
|||
m_T.topRows(i+2).applyOnTheRight(i+1,i,G); |
|||
// update Z |
|||
if (m_computeQZ) |
|||
m_Z.applyOnTheLeft(i+1,i,G.adjoint()); |
|||
|
|||
m_S.coeffRef(i+1,i) = Scalar(0.0); |
|||
m_T.coeffRef(i+1,i) = Scalar(0.0); |
|||
} |
|||
} |
|||
else |
|||
{ |
|||
pushDownZero(z,i,i+1); |
|||
} |
|||
} |
|||
|
|||
/** \internal use zero in T(z,z) to zero S(l,l-1), working in block f..l */ |
|||
template<typename MatrixType> |
|||
inline void RealQZ<MatrixType>::pushDownZero(Index z, Index f, Index l) |
|||
{ |
|||
JRs G; |
|||
const Index dim = m_S.cols(); |
|||
for (Index zz=z; zz<l; zz++) |
|||
{ |
|||
// push 0 down |
|||
Index firstColS = zz>f ? (zz-1) : zz; |
|||
G.makeGivens(m_T.coeff(zz, zz+1), m_T.coeff(zz+1, zz+1)); |
|||
m_S.rightCols(dim-firstColS).applyOnTheLeft(zz,zz+1,G.adjoint()); |
|||
m_T.rightCols(dim-zz).applyOnTheLeft(zz,zz+1,G.adjoint()); |
|||
m_T.coeffRef(zz+1,zz+1) = Scalar(0.0); |
|||
// update Q |
|||
if (m_computeQZ) |
|||
m_Q.applyOnTheRight(zz,zz+1,G); |
|||
// kill S(zz+1, zz-1) |
|||
if (zz>f) |
|||
{ |
|||
G.makeGivens(m_S.coeff(zz+1, zz), m_S.coeff(zz+1,zz-1)); |
|||
m_S.topRows(zz+2).applyOnTheRight(zz, zz-1,G); |
|||
m_T.topRows(zz+1).applyOnTheRight(zz, zz-1,G); |
|||
m_S.coeffRef(zz+1,zz-1) = Scalar(0.0); |
|||
// update Z |
|||
if (m_computeQZ) |
|||
m_Z.applyOnTheLeft(zz,zz-1,G.adjoint()); |
|||
} |
|||
} |
|||
// finally kill S(l,l-1) |
|||
G.makeGivens(m_S.coeff(l,l), m_S.coeff(l,l-1)); |
|||
m_S.applyOnTheRight(l,l-1,G); |
|||
m_T.applyOnTheRight(l,l-1,G); |
|||
m_S.coeffRef(l,l-1)=Scalar(0.0); |
|||
// update Z |
|||
if (m_computeQZ) |
|||
m_Z.applyOnTheLeft(l,l-1,G.adjoint()); |
|||
} |
|||
|
|||
/** \internal QR-like iterative step for block f..l */ |
|||
template<typename MatrixType> |
|||
inline void RealQZ<MatrixType>::step(Index f, Index l, Index iter) { |
|||
const Index dim = m_S.cols(); |
|||
|
|||
// x, y, z |
|||
Scalar x, y, z; |
|||
if (iter==10) |
|||
{ |
|||
// Wilkinson ad hoc shift |
|||
const Scalar |
|||
a11=m_S.coeff(f+0,f+0), a12=m_S.coeff(f+0,f+1), |
|||
a21=m_S.coeff(f+1,f+0), a22=m_S.coeff(f+1,f+1), a32=m_S.coeff(f+2,f+1), |
|||
b12=m_T.coeff(f+0,f+1), |
|||
b11i=Scalar(1.0)/m_T.coeff(f+0,f+0), |
|||
b22i=Scalar(1.0)/m_T.coeff(f+1,f+1), |
|||
a87=m_S.coeff(l-1,l-2), |
|||
a98=m_S.coeff(l-0,l-1), |
|||
b77i=Scalar(1.0)/m_T.coeff(l-2,l-2), |
|||
b88i=Scalar(1.0)/m_T.coeff(l-1,l-1); |
|||
Scalar ss = internal::abs(a87*b77i) + internal::abs(a98*b88i), |
|||
lpl = Scalar(1.5)*ss, |
|||
ll = ss*ss; |
|||
x = ll + a11*a11*b11i*b11i - lpl*a11*b11i + a12*a21*b11i*b22i |
|||
- a11*a21*b12*b11i*b11i*b22i; |
|||
y = a11*a21*b11i*b11i - lpl*a21*b11i + a21*a22*b11i*b22i |
|||
- a21*a21*b12*b11i*b11i*b22i; |
|||
z = a21*a32*b11i*b22i; |
|||
} |
|||
else if (iter==16) |
|||
{ |
|||
// another exceptional shift |
|||
x = m_S.coeff(f,f)/m_T.coeff(f,f)-m_S.coeff(l,l)/m_T.coeff(l,l) + m_S.coeff(l,l-1)*m_T.coeff(l-1,l) / |
|||
(m_T.coeff(l-1,l-1)*m_T.coeff(l,l)); |
|||
y = m_S.coeff(f+1,f)/m_T.coeff(f,f); |
|||
z = 0; |
|||
} |
|||
else if (iter>23 && !(iter%8)) |
|||
{ |
|||
// extremely exceptional shift |
|||
x = internal::random<Scalar>(-1.0,1.0); |
|||
y = internal::random<Scalar>(-1.0,1.0); |
|||
z = internal::random<Scalar>(-1.0,1.0); |
|||
} |
|||
else |
|||
{ |
|||
// Compute the shifts: (x,y,z,0...) = (AB^-1 - l1 I) (AB^-1 - l2 I) e1 |
|||
// where l1 and l2 are the eigenvalues of the 2x2 matrix C = U V^-1 where |
|||
// U and V are 2x2 bottom right sub matrices of A and B. Thus: |
|||
// = AB^-1AB^-1 + l1 l2 I - (l1+l2)(AB^-1) |
|||
// = AB^-1AB^-1 + det(M) - tr(M)(AB^-1) |
|||
// Since we are only interested in having x, y, z with a correct ratio, we have: |
|||
const Scalar |
|||
a11 = m_S.coeff(f,f), a12 = m_S.coeff(f,f+1), |
|||
a21 = m_S.coeff(f+1,f), a22 = m_S.coeff(f+1,f+1), |
|||
a32 = m_S.coeff(f+2,f+1), |
|||
|
|||
a88 = m_S.coeff(l-1,l-1), a89 = m_S.coeff(l-1,l), |
|||
a98 = m_S.coeff(l,l-1), a99 = m_S.coeff(l,l), |
|||
|
|||
b11 = m_T.coeff(f,f), b12 = m_T.coeff(f,f+1), |
|||
b22 = m_T.coeff(f+1,f+1), |
|||
|
|||
b88 = m_T.coeff(l-1,l-1), b89 = m_T.coeff(l-1,l), |
|||
b99 = m_T.coeff(l,l); |
|||
|
|||
x = ( (a88/b88 - a11/b11)*(a99/b99 - a11/b11) - (a89/b99)*(a98/b88) + (a98/b88)*(b89/b99)*(a11/b11) ) * (b11/a21) |
|||
+ a12/b22 - (a11/b11)*(b12/b22); |
|||
y = (a22/b22-a11/b11) - (a21/b11)*(b12/b22) - (a88/b88-a11/b11) - (a99/b99-a11/b11) + (a98/b88)*(b89/b99); |
|||
z = a32/b22; |
|||
} |
|||
|
|||
JRs G; |
|||
|
|||
for (Index k=f; k<=l-2; k++) |
|||
{ |
|||
// variables for Householder reflections |
|||
Vector2s essential2; |
|||
Scalar tau, beta; |
|||
|
|||
Vector3s hr(x,y,z); |
|||
|
|||
// Q_k to annihilate S(k+1,k-1) and S(k+2,k-1) |
|||
hr.makeHouseholderInPlace(tau, beta); |
|||
essential2 = hr.template bottomRows<2>(); |
|||
Index fc=(std::max)(k-1,Index(0)); // first col to update |
|||
m_S.template middleRows<3>(k).rightCols(dim-fc).applyHouseholderOnTheLeft(essential2, tau, m_workspace.data()); |
|||
m_T.template middleRows<3>(k).rightCols(dim-fc).applyHouseholderOnTheLeft(essential2, tau, m_workspace.data()); |
|||
if (m_computeQZ) |
|||
m_Q.template middleCols<3>(k).applyHouseholderOnTheRight(essential2, tau, m_workspace.data()); |
|||
if (k>f) |
|||
m_S.coeffRef(k+2,k-1) = m_S.coeffRef(k+1,k-1) = Scalar(0.0); |
|||
|
|||
// Z_{k1} to annihilate T(k+2,k+1) and T(k+2,k) |
|||
hr << m_T.coeff(k+2,k+2),m_T.coeff(k+2,k),m_T.coeff(k+2,k+1); |
|||
hr.makeHouseholderInPlace(tau, beta); |
|||
essential2 = hr.template bottomRows<2>(); |
|||
{ |
|||
Index lr = (std::min)(k+4,dim); // last row to update |
|||
Map<Matrix<Scalar,Dynamic,1> > tmp(m_workspace.data(),lr); |
|||
// S |
|||
tmp = m_S.template middleCols<2>(k).topRows(lr) * essential2; |
|||
tmp += m_S.col(k+2).head(lr); |
|||
m_S.col(k+2).head(lr) -= tau*tmp; |
|||
m_S.template middleCols<2>(k).topRows(lr) -= (tau*tmp) * essential2.adjoint(); |
|||
// T |
|||
tmp = m_T.template middleCols<2>(k).topRows(lr) * essential2; |
|||
tmp += m_T.col(k+2).head(lr); |
|||
m_T.col(k+2).head(lr) -= tau*tmp; |
|||
m_T.template middleCols<2>(k).topRows(lr) -= (tau*tmp) * essential2.adjoint(); |
|||
} |
|||
if (m_computeQZ) |
|||
{ |
|||
// Z |
|||
Map<Matrix<Scalar,1,Dynamic> > tmp(m_workspace.data(),dim); |
|||
tmp = essential2.adjoint()*(m_Z.template middleRows<2>(k)); |
|||
tmp += m_Z.row(k+2); |
|||
m_Z.row(k+2) -= tau*tmp; |
|||
m_Z.template middleRows<2>(k) -= essential2 * (tau*tmp); |
|||
} |
|||
m_T.coeffRef(k+2,k) = m_T.coeffRef(k+2,k+1) = Scalar(0.0); |
|||
|
|||
// Z_{k2} to annihilate T(k+1,k) |
|||
G.makeGivens(m_T.coeff(k+1,k+1), m_T.coeff(k+1,k)); |
|||
m_S.applyOnTheRight(k+1,k,G); |
|||
m_T.applyOnTheRight(k+1,k,G); |
|||
// update Z |
|||
if (m_computeQZ) |
|||
m_Z.applyOnTheLeft(k+1,k,G.adjoint()); |
|||
m_T.coeffRef(k+1,k) = Scalar(0.0); |
|||
|
|||
// update x,y,z |
|||
x = m_S.coeff(k+1,k); |
|||
y = m_S.coeff(k+2,k); |
|||
if (k < l-2) |
|||
z = m_S.coeff(k+3,k); |
|||
} // loop over k |
|||
|
|||
// Q_{n-1} to annihilate y = S(l,l-2) |
|||
G.makeGivens(x,y); |
|||
m_S.applyOnTheLeft(l-1,l,G.adjoint()); |
|||
m_T.applyOnTheLeft(l-1,l,G.adjoint()); |
|||
if (m_computeQZ) |
|||
m_Q.applyOnTheRight(l-1,l,G); |
|||
m_S.coeffRef(l,l-2) = Scalar(0.0); |
|||
|
|||
// Z_{n-1} to annihilate T(l,l-1) |
|||
G.makeGivens(m_T.coeff(l,l),m_T.coeff(l,l-1)); |
|||
m_S.applyOnTheRight(l,l-1,G); |
|||
m_T.applyOnTheRight(l,l-1,G); |
|||
if (m_computeQZ) |
|||
m_Z.applyOnTheLeft(l,l-1,G.adjoint()); |
|||
m_T.coeffRef(l,l-1) = Scalar(0.0); |
|||
} |
|||
|
|||
|
|||
template<typename MatrixType> |
|||
RealQZ<MatrixType>& RealQZ<MatrixType>::compute(const MatrixType& A_in, const MatrixType& B_in, bool computeQZ) |
|||
{ |
|||
|
|||
const Index dim = A_in.cols(); |
|||
|
|||
assert (A_in.rows()==dim && A_in.cols()==dim |
|||
&& B_in.rows()==dim && B_in.cols()==dim |
|||
&& "Need square matrices of the same dimension"); |
|||
|
|||
m_isInitialized = true; |
|||
m_computeQZ = computeQZ; |
|||
m_S = A_in; m_T = B_in; |
|||
m_workspace.resize(dim*2); |
|||
m_global_iter = 0; |
|||
|
|||
// entrance point: hessenberg triangular decomposition |
|||
hessenbergTriangular(); |
|||
// compute L1 vector norms of T, S into m_normOfS, m_normOfT |
|||
computeNorms(); |
|||
|
|||
Index l = dim-1, |
|||
f, |
|||
local_iter = 0; |
|||
|
|||
while (l>0 && local_iter<m_maxIters) |
|||
{ |
|||
f = findSmallSubdiagEntry(l); |
|||
// now rows and columns f..l (including) decouple from the rest of the problem |
|||
if (f>0) m_S.coeffRef(f,f-1) = Scalar(0.0); |
|||
if (f == l) // One root found |
|||
{ |
|||
l--; |
|||
local_iter = 0; |
|||
} |
|||
else if (f == l-1) // Two roots found |
|||
{ |
|||
splitOffTwoRows(f); |
|||
l -= 2; |
|||
local_iter = 0; |
|||
} |
|||
else // No convergence yet |
|||
{ |
|||
// if there's zero on diagonal of T, we can isolate an eigenvalue with Givens rotations |
|||
Index z = findSmallDiagEntry(f,l); |
|||
if (z>=f) |
|||
{ |
|||
// zero found |
|||
pushDownZero(z,f,l); |
|||
} |
|||
else |
|||
{ |
|||
// We are sure now that S.block(f,f, l-f+1,l-f+1) is underuced upper-Hessenberg |
|||
// and T.block(f,f, l-f+1,l-f+1) is invertible uper-triangular, which allows to |
|||
// apply a QR-like iteration to rows and columns f..l. |
|||
step(f,l, local_iter); |
|||
local_iter++; |
|||
m_global_iter++; |
|||
} |
|||
} |
|||
} |
|||
// check if we converged before reaching iterations limit |
|||
m_info = (local_iter<m_maxIters) ? Success : NoConvergence; |
|||
return *this; |
|||
} // end compute |
|||
|
|||
} // end namespace Eigen |
|||
|
|||
#endif //EIGEN_REAL_QZ |
@ -0,0 +1,6 @@ |
|||
FILE(GLOB Eigen_MetisSupport_SRCS "*.h") |
|||
|
|||
INSTALL(FILES |
|||
${Eigen_MetisSupport_SRCS} |
|||
DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/MetisSupport COMPONENT Devel |
|||
) |
@ -0,0 +1,138 @@ |
|||
// This file is part of Eigen, a lightweight C++ template library |
|||
// for linear algebra. |
|||
// |
|||
// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr> |
|||
// |
|||
// This Source Code Form is subject to the terms of the Mozilla |
|||
// Public License v. 2.0. If a copy of the MPL was not distributed |
|||
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. |
|||
#ifndef METIS_SUPPORT_H |
|||
#define METIS_SUPPORT_H |
|||
|
|||
namespace Eigen { |
|||
/** |
|||
* Get the fill-reducing ordering from the METIS package |
|||
* |
|||
* If A is the original matrix and Ap is the permuted matrix, |
|||
* the fill-reducing permutation is defined as follows : |
|||
* Row (column) i of A is the matperm(i) row (column) of Ap. |
|||
* WARNING: As computed by METIS, this corresponds to the vector iperm (instead of perm) |
|||
*/ |
|||
template <typename Index> |
|||
class MetisOrdering |
|||
{ |
|||
public: |
|||
typedef PermutationMatrix<Dynamic,Dynamic,Index> PermutationType; |
|||
typedef Matrix<Index,Dynamic,1> IndexVector; |
|||
|
|||
template <typename MatrixType> |
|||
void get_symmetrized_graph(const MatrixType& A) |
|||
{ |
|||
Index m = A.cols(); |
|||
|
|||
// Get the transpose of the input matrix |
|||
MatrixType At = A.transpose(); |
|||
// Get the number of nonzeros elements in each row/col of At+A |
|||
Index TotNz = 0; |
|||
IndexVector visited(m); |
|||
visited.setConstant(-1); |
|||
for (int j = 0; j < m; j++) |
|||
{ |
|||
// Compute the union structure of of A(j,:) and At(j,:) |
|||
visited(j) = j; // Do not include the diagonal element |
|||
// Get the nonzeros in row/column j of A |
|||
for (typename MatrixType::InnerIterator it(A, j); it; ++it) |
|||
{ |
|||
Index idx = it.index(); // Get the row index (for column major) or column index (for row major) |
|||
if (visited(idx) != j ) |
|||
{ |
|||
visited(idx) = j; |
|||
++TotNz; |
|||
} |
|||
} |
|||
//Get the nonzeros in row/column j of At |
|||
for (typename MatrixType::InnerIterator it(At, j); it; ++it) |
|||
{ |
|||
Index idx = it.index(); |
|||
if(visited(idx) != j) |
|||
{ |
|||
visited(idx) = j; |
|||
++TotNz; |
|||
} |
|||
} |
|||
} |
|||
// Reserve place for A + At |
|||
m_indexPtr.resize(m+1); |
|||
m_innerIndices.resize(TotNz); |
|||
|
|||
// Now compute the real adjacency list of each column/row |
|||
visited.setConstant(-1); |
|||
Index CurNz = 0; |
|||
for (int j = 0; j < m; j++) |
|||
{ |
|||
m_indexPtr(j) = CurNz; |
|||
|
|||
visited(j) = j; // Do not include the diagonal element |
|||
// Add the pattern of row/column j of A to A+At |
|||
for (typename MatrixType::InnerIterator it(A,j); it; ++it) |
|||
{ |
|||
Index idx = it.index(); // Get the row index (for column major) or column index (for row major) |
|||
if (visited(idx) != j ) |
|||
{ |
|||
visited(idx) = j; |
|||
m_innerIndices(CurNz) = idx; |
|||
CurNz++; |
|||
} |
|||
} |
|||
//Add the pattern of row/column j of At to A+At |
|||
for (typename MatrixType::InnerIterator it(At, j); it; ++it) |
|||
{ |
|||
Index idx = it.index(); |
|||
if(visited(idx) != j) |
|||
{ |
|||
visited(idx) = j; |
|||
m_innerIndices(CurNz) = idx; |
|||
++CurNz; |
|||
} |
|||
} |
|||
} |
|||
m_indexPtr(m) = CurNz; |
|||
} |
|||
|
|||
template <typename MatrixType> |
|||
void operator() (const MatrixType& A, PermutationType& matperm) |
|||
{ |
|||
Index m = A.cols(); |
|||
IndexVector perm(m),iperm(m); |
|||
// First, symmetrize the matrix graph. |
|||
get_symmetrized_graph(A); |
|||
int output_error; |
|||
|
|||
// Call the fill-reducing routine from METIS |
|||
output_error = METIS_NodeND(&m, m_indexPtr.data(), m_innerIndices.data(), NULL, NULL, perm.data(), iperm.data()); |
|||
|
|||
if(output_error != METIS_OK) |
|||
{ |
|||
//FIXME The ordering interface should define a class of possible errors |
|||
std::cerr << "ERROR WHILE CALLING THE METIS PACKAGE \n"; |
|||
return; |
|||
} |
|||
|
|||
// Get the fill-reducing permutation |
|||
//NOTE: If Ap is the permuted matrix then perm and iperm vectors are defined as follows |
|||
// Row (column) i of Ap is the perm(i) row(column) of A, and row (column) i of A is the iperm(i) row(column) of Ap |
|||
|
|||
// To be consistent with the use of the permutation in SparseLU module, we thus keep the iperm vector |
|||
matperm.resize(m); |
|||
for (int j = 0; j < m; j++) |
|||
matperm.indices()(j) = iperm(j); |
|||
|
|||
} |
|||
|
|||
protected: |
|||
IndexVector m_indexPtr; // Pointer to the adjacenccy list of each row/column |
|||
IndexVector m_innerIndices; // Adjacency list |
|||
}; |
|||
|
|||
}// end namespace eigen |
|||
#endif |
1849
resources/3rdparty/eigen/Eigen/src/OrderingMethods/Eigen_Colamd.h
File diff suppressed because it is too large
View File
File diff suppressed because it is too large
View File
@ -0,0 +1,158 @@ |
|||
|
|||
// This file is part of Eigen, a lightweight C++ template library |
|||
// for linear algebra. |
|||
// |
|||
// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr> |
|||
// |
|||
// Eigen is free software; you can redistribute it and/or |
|||
// modify it under the terms of the GNU Lesser General Public |
|||
// License as published by the Free Software Foundation; either |
|||
// version 3 of the License, or (at your option) any later version. |
|||
// |
|||
// Alternatively, you can redistribute it and/or |
|||
// modify it under the terms of the GNU General Public License as |
|||
// published by the Free Software Foundation; either version 2 of |
|||
// the License, or (at your option) any later version. |
|||
// |
|||
// Eigen 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 Lesser General Public License or the |
|||
// GNU General Public License for more details. |
|||
// |
|||
// You should have received a copy of the GNU Lesser General Public |
|||
// License and a copy of the GNU General Public License along with |
|||
// Eigen. If not, see <http://www.gnu.org/licenses/>. |
|||
|
|||
#ifndef EIGEN_ORDERING_H |
|||
#define EIGEN_ORDERING_H |
|||
|
|||
#include "Amd.h" |
|||
namespace Eigen { |
|||
|
|||
#include "Eigen_Colamd.h" |
|||
|
|||
namespace internal { |
|||
|
|||
/** |
|||
* Get the symmetric pattern A^T+A from the input matrix A. |
|||
* FIXME: The values should not be considered here |
|||
*/ |
|||
template<typename MatrixType> |
|||
void ordering_helper_at_plus_a(const MatrixType& mat, MatrixType& symmat) |
|||
{ |
|||
MatrixType C; |
|||
C = mat.transpose(); // NOTE: Could be costly |
|||
for (int i = 0; i < C.rows(); i++) |
|||
{ |
|||
for (typename MatrixType::InnerIterator it(C, i); it; ++it) |
|||
it.valueRef() = 0.0; |
|||
} |
|||
symmat = C + mat; |
|||
} |
|||
|
|||
} |
|||
|
|||
/** |
|||
* Get the approximate minimum degree ordering |
|||
* If the matrix is not structurally symmetric, an ordering of A^T+A is computed |
|||
* \tparam Index The type of indices of the matrix |
|||
*/ |
|||
template <typename Index> |
|||
class AMDOrdering |
|||
{ |
|||
public: |
|||
typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType; |
|||
|
|||
/** Compute the permutation vector from a sparse matrix |
|||
* This routine is much faster if the input matrix is column-major |
|||
*/ |
|||
template <typename MatrixType> |
|||
void operator()(const MatrixType& mat, PermutationType& perm) |
|||
{ |
|||
// Compute the symmetric pattern |
|||
SparseMatrix<typename MatrixType::Scalar, ColMajor, Index> symm; |
|||
internal::ordering_helper_at_plus_a(mat,symm); |
|||
|
|||
// Call the AMD routine |
|||
//m_mat.prune(keep_diag()); |
|||
internal::minimum_degree_ordering(symm, perm); |
|||
} |
|||
|
|||
/** Compute the permutation with a selfadjoint matrix */ |
|||
template <typename SrcType, unsigned int SrcUpLo> |
|||
void operator()(const SparseSelfAdjointView<SrcType, SrcUpLo>& mat, PermutationType& perm) |
|||
{ |
|||
SparseMatrix<typename SrcType::Scalar, ColMajor, Index> C = mat; |
|||
|
|||
// Call the AMD routine |
|||
// m_mat.prune(keep_diag()); //Remove the diagonal elements |
|||
internal::minimum_degree_ordering(C, perm); |
|||
} |
|||
}; |
|||
|
|||
/** |
|||
* Get the natural ordering |
|||
* |
|||
*NOTE Returns an empty permutation matrix |
|||
* \tparam Index The type of indices of the matrix |
|||
*/ |
|||
template <typename Index> |
|||
class NaturalOrdering |
|||
{ |
|||
public: |
|||
typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType; |
|||
|
|||
/** Compute the permutation vector from a column-major sparse matrix */ |
|||
template <typename MatrixType> |
|||
void operator()(const MatrixType& mat, PermutationType& perm) |
|||
{ |
|||
perm.resize(0); |
|||
} |
|||
|
|||
}; |
|||
|
|||
/** |
|||
* Get the column approximate minimum degree ordering |
|||
* The matrix should be in column-major format |
|||
*/ |
|||
template<typename Index> |
|||
class COLAMDOrdering; |
|||
#include "Eigen_Colamd.h" |
|||
|
|||
template<typename Index> |
|||
class COLAMDOrdering |
|||
{ |
|||
public: |
|||
typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType; |
|||
typedef Matrix<Index, Dynamic, 1> IndexVector; |
|||
/** Compute the permutation vector form a sparse matrix */ |
|||
template <typename MatrixType> |
|||
void operator() (const MatrixType& mat, PermutationType& perm) |
|||
{ |
|||
int m = mat.rows(); |
|||
int n = mat.cols(); |
|||
int nnz = mat.nonZeros(); |
|||
// Get the recommended value of Alen to be used by colamd |
|||
int Alen = internal::colamd_recommended(nnz, m, n); |
|||
// Set the default parameters |
|||
double knobs [COLAMD_KNOBS]; |
|||
int stats [COLAMD_STATS]; |
|||
internal::colamd_set_defaults(knobs); |
|||
|
|||
int info; |
|||
IndexVector p(n+1), A(Alen); |
|||
for(int i=0; i <= n; i++) p(i) = mat.outerIndexPtr()[i]; |
|||
for(int i=0; i < nnz; i++) A(i) = mat.innerIndexPtr()[i]; |
|||
// Call Colamd routine to compute the ordering |
|||
info = internal::colamd(m, n, Alen, A.data(), p.data(), knobs, stats); |
|||
eigen_assert( info && "COLAMD failed " ); |
|||
|
|||
perm.resize(n); |
|||
for (int i = 0; i < n; i++) perm.indices()(p(i)) = i; |
|||
|
|||
} |
|||
|
|||
}; |
|||
|
|||
} // end namespace Eigen |
|||
#endif |
Some files were not shown because too many files changed in this diff
Reference in new issue
xxxxxxxxxx