Log of /trunk/math/matrix/inc/TDecompLU.h
Parent Directory
Revision
22885 -
(
view)
(
download)
(
as text)
(
annotate)
-
[select for diffs]
Modified
Fri Mar 28 13:57:25 2008 UTC (6 years, 9 months ago) by
rdm
File length: 3570 byte(s)
Diff to
previous 20882
move the following directories under the new "math" meta directory:
mathcore
mathmore
fftw
foam
fumili
genvector
matrix
minuit
minuit2
mlp
physics
smatrix
splot
unuran
quadp
Revision
16458 -
(
view)
(
download)
(
as text)
(
annotate)
-
[select for diffs]
Modified
Fri Oct 6 06:52:34 2006 UTC (8 years, 3 months ago) by
brun
Original Path:
trunk/matrix/inc/TDecompLU.h
File length: 3630 byte(s)
Diff to
previous 15174
From Eddy Offermann:
Previously, the matrix package contained a large amount of ASSERT
statements which are a nuisance for programs analyzing a series of
independent events like in high-energy physics .
Assert's were issued when for instance a division by zero was requested
or a matrix was invalid . Most algorithms made matrices/vectors invalid
after an error occured in an operation , like inversion of a singular matrix .
Unfortunately, not all assert's were accompanied by error messages .
This situation has been completely overhauled :
- All error conditions in the algorithms are now accompanied by error
messages .
- In all algorithms it is still asserted that vectors/matrices are
valid BUT only in very few cases is a matrix/vector made invalid :
for instance if memory is allocated with incorrect parameters .
- In case of division by zero, the division is skipped . In case of a
singular matrix, the inversion routine returns the original matrix .
In the past the result of an inversion could be checked through the
value of the returned determinant or checking whether the inverted matrix
was valid .
Since from now on, we never make the matrix invalid in this operation, the
latter check will not indicate a singular matrix anymore .
The decompostion classes TDecomp... have a backward-compatible
change in the interface which makes detection of singularity easier :
old interface :
void Invert (TMatrixD &inv);
TMatrixD Invert ();
new interface :
Bool_t Invert (TMatrixD &inv);
TMatrixD Invert (Bool_t &status);
TMatrixD Invert () { Bool_t status; return Invert(status);
}
The returned status is kFALSE in case of singularity .
The old situation is easily reproduced by setting the ROOT variable
gErrorAbortLevel to kError . This cause an exception when there is an error
message (Error...) ., In the past the matrix would be made invalid which
would cause the next operation to throw an exception .
Revision
10329 -
(
view)
(
download)
(
as text)
(
annotate)
-
[select for diffs]
Modified
Sat Oct 16 18:09:17 2004 UTC (10 years, 3 months ago) by
brun
Original Path:
trunk/matrix/inc/TDecompLU.h
File length: 3505 byte(s)
Diff to
previous 9503
From Eddy Offermann:
- added new decomposition class TDecompBK: Bunch-Kaufman algorithm
this class is designed to decompose real symmetric matrices . Now
we can guarantee that inverse of a symmetric matrix is again
symmetric .
- As a consequence of the introduction of TDecompBK, we added
Invert(..) functionality to TMatrixD/FSym
Also added InvertFast to TMatrixD/FSym which uses Cramer
inversion for matrices <= (6 x 6) , resulting in large speedup
but less accuracy . The necessary routines are encapsulated
in the new TMatrixDSymCramerInv class .
Note of caution :
Unlike the other decompositions, Bunch-Kaufmann does not result
in a triangular matrix of which the determinant is trivially calculated
by multiplying the diagonal elements . As a consequence,
we have to invoke yet another decomposition when the user
requests the determinant .
So Invert() will be roughly twice faster than Double_t *det; Invert(det)
- Changed the algoritm for the Cholesky decomposition so that we
do not need to have a copy of the original matrix around !
After removal of class member fA, increased the version to 2
- Added sorting of eigen -values/-vectors to TMatrixDEigen .
They are now sorted according to Re^2+Im^2 of the (possible)
complex eigenvalues in descending order . TMatrixDSymEigen
already sorted . Since most users use TMatrixD instead of
TMatrixDSym, this will ensure backward compatiblity since
old root contained the symmetric algorithms in TMatrixD .
- Added rank 1 update algorihms : A += alpha * x y^T
(x and y are vectors)
to the non-sparse matrix classes (adding this to sparse
would make this matrix immediately dense :)
- Moved Invert fom TDecompBase to the different decomposition
classes . Reason is that some classes can decompose A(m,n)
where m > n . The result of Invert() is then a "pseudo-invert"
: A_pseudo = (A^T A)^-1 A^T .
- Removed an unfortunate bug in the TMatrixD/FSub class :
both operators
void operator+=(Double_t val);
void operator*=(Double_t val);
performed actually
void operator=(Double_t val);
~
Revision
9503 -
(
view)
(
download)
(
as text)
(
annotate)
-
[select for diffs]
Modified
Mon Jul 12 20:00:41 2004 UTC (10 years, 6 months ago) by
brun
Original Path:
trunk/matrix/inc/TDecompLU.h
File length: 3366 byte(s)
Diff to
previous 9181
From Eddy Offermann:
Some code re-organization:
- made TDecompLU::DecomposeLUCrout and TDecompLU::DecomposeLUGauss
protected.
As a consequence changed TDecompLU::InvertLU so that it now calls
the decomposition step .
This change prohibits the user of calling InvertLU without applying
first
the decomposition step
- Routines in TMatrixD/FCramerInv return now a Bool_t instead of Int_t
.
More important, these routines check that the matrix is indeed of the
appropriate shape, since the user could call them directly instead
of going through TMatrixD/F::InvertFast
- Removed the redundant TDecompBase::MultiSolve(TMatrixDSym &B), it did
exactly what TDecompBase::MultiSolve(TMatrixD &B) did .
- void (Trans)Solve(TVectorD &b,Int_t ok) and
void (Trans)Solve(TMatrixDColumn &cb, Int_t ok) routines in TDecomp*
now make vector b or matrix b in case of an error (when ok was set to
false)
Revision
9032 -
(
view)
(
download)
(
as text)
(
annotate)
-
[select for diffs]
Modified
Thu May 27 06:39:53 2004 UTC (10 years, 8 months ago) by
brun
Original Path:
trunk/matrix/inc/TDecompLU.h
File length: 3425 byte(s)
Diff to
previous 8892
From Eddy Offermann
Implemented the matrix status through TObject::fBits .
Also removed the fStatus data member of TDecompbase
and use now TObject::fBits .
User can now reset the matrix status through
MakeValid()
Revision
8892 -
(
view)
(
download)
(
as text)
(
annotate)
-
[select for diffs]
Modified
Wed May 12 10:39:29 2004 UTC (10 years, 8 months ago) by
brun
Original Path:
trunk/matrix/inc/TDecompLU.h
File length: 3430 byte(s)
Diff to
previous 8456
From Eddy Offermann:
1) patch for the matrix directory (cvs diff -a matrix)
2) tar file containing new modules:
matrix/inc/TDecompSparse.h
matrix/inc/TMatrixDSparse.h
matrix/src/TDecompSparse.cxx
matrix/src/TMatrixDSparse.cxx
3) new stressLinear.cxx
Description of changes:
1) many updates of comments in code
2) possibility to construct random matrices, general
symmetric and positive definite
3) Introduction of new matrix class: TMatrixDSparse .
It defines a general sparse matrix in
Harwell-Boeing
format (sparse structure definition optimized for
matrix-vector multiplications).
Like all the other matrix classes, it derives from
TMatrixDBase
4) Introduction of a new decomposition class:
TDecompSparse . It allows decomposition (and
therefore equation solving) of sparse symmetric
matrices
Revision
8456 -
(
view)
(
download)
(
as text)
(
annotate)
-
[select for diffs]
Modified
Mon Mar 22 08:34:36 2004 UTC (10 years, 10 months ago) by
brun
Original Path:
trunk/matrix/inc/TDecompLU.h
File length: 3426 byte(s)
Diff to
previous 8442
From Eddy:
Added to the decomposition class constructors in
which just the matrix size is given as argument .
Also a new state (kMatrixSet) was added .
Instead of
TMatrixD a;
TDecompLU lu(a);
one can do
TMatrixD a;
TDecompLU lu(n);
lu.SetMatrix(a);
The SetMatrix command resizes, if necessary, the
data members in the decomposition class .
This set up willl be usefull when one wants to
decompose
matrices in a loop
Revision
8442 -
(
view)
(
download)
(
as text)
(
annotate)
-
[select for diffs]
Modified
Fri Mar 19 14:20:40 2004 UTC (10 years, 10 months ago) by
brun
Original Path:
trunk/matrix/inc/TDecompLU.h
File length: 3349 byte(s)
Diff to
previous 8225
From Eddy Offermann:
1. Removed streaming of matrix view classes (TMatrixDRow,TMatrixDColumn,TMatrixDDiag,TMatrixDFlat)
because a. we want these classes to be as light as possible
b. old implementation was wrong, some member elements were "(!)" out
2. Added a series of functions to support upcoming quadratic programming functionality:
TVectorD/F:
- AddElemMult(TVectorF &t,Float_t alpha,const TVectorF &s1,const TVectorF &s2);
t += alpha * s1*s2 , where "s1*s2 is an element-wise multiplication
- AddElemMult(TVectorF &t,Float_t alpha,const TVectorF &s1,const TVectorF &s2,const TVectorF &map);
t += alpha * s1*s2 , where "s1*s2 is an element-wise multiplication. It is only performed
for those elements i where map[i] != 0
- AddElemDiv (TVectorF &t,Float_t alpha,const TVectorF &s1,const TVectorF &s2);
t += alpha * s1\s2 , where "s1\s2 is an element-wise division
- AddElemDiv (TVectorF &t,Float_t alpha,const TVectorF &s1,const TVectorF &s2,const TVectorF &map);
t += alpha * s1\s2 , where "s1\s2 is an element-wise division. It is only performed
for those elements i where map[i] != 0
- ElementMult(TVectorF &t,const TVectorF &s,const TVectorF &map);
t += t*s , where "t*s" is an element-wise multiplication. It is only performed
for those elements i where map[i] != 0
- ElementDiv (TVectorF &t,const TVectorF &s,const TVectorF &map);
t += t/s , where "t/s" is an element-wise division. It is only performed
for those elements i where map[i] != 0
- TVectorD &SelectNonZeros (const TVectorD &select);
put elements for which select[i] == 0 to zero
- Bool_t MatchesNonZeroPattern(const TVectorD &select);
verify that elements are != 0 for which select[i] != 0
- Bool_t SomePositive (const TVectorD &select);
verify that all elements > 0 for which select[i] != 0
- void AddSomeConstant (Double_t val,const TVectorD &select);
Add val to those elements for which select[i] != 0
- TVectorF &TVectorF::Invert()
v[i] = 1/v[i]
TVectorD/F and TMatrixD/FBase:
- virtual Double_t Sum () const;
- virtual Double_t Min () const;
- virtual Double_t Max () const;
3. Added another data member fA to TDecompChol . Unlike the other decomposition schemes, it is not
possible to store the original input matrix in one of the decomposition members
4. Added method "SetMatrix(const TMatrix.. &)" to all decomposition classes . It allows to
change the matrix that should be decomposed. Of course it resets the class status to kInit
and triggers a new factorization.
5. Added private methods to TMatrixD and TMatrixDSym
const TMatrixD EigenVectors(TVectorD &eigenValues) const; // This function is now obsolete (and not implemented), you
// should use TMatrixDSymEigen or TMatrixDEigen .
6. Changed the name "Adopt" to "Use" to be in agreement with the convention that the memory management
of adopted members is now the reponsibility of the class
Revision
8225 -
(
view)
(
download)
(
as text)
(
annotate)
-
[select for diffs]
Modified
Wed Feb 18 14:36:43 2004 UTC (10 years, 11 months ago) by
brun
Original Path:
trunk/matrix/inc/TDecompLU.h
File length: 3291 byte(s)
Diff to
previous 8189
From Eddy Offermann:
First, a patch to the matrix directory is attached. It
turned out that that 3x3 Jacobi inversion code, did
not
calculate the determinant right . Of course,
stressLinear caught and I did not run it :)
Revision
8118 -
(
view)
(
download)
(
as text)
(
annotate)
-
[select for diffs]
Modified
Wed Feb 4 17:12:44 2004 UTC (10 years, 11 months ago) by
brun
Original Path:
trunk/matrix/inc/TDecompLU.h
File length: 2883 byte(s)
Diff to
previous 8112
From Eddy Offermann:
1. LinkDef.h : Added NormalEqn friend functions from
TDecompChol class
2. TMatrixDBase: GetDecompMatrix() returns a TMatrixD
instead of TMatrixDBase
Added Solve and TransSolve functions
that return the solution (others
replace
the input argument which is
inconvinient
when output vector size is smaller
than
input (in ncase of fit).
Also added Invert that returns matrix
instead of replacing the input
argument.
This is convinient for a generalized
(pseudo) inverse: in (m x n) out (n x
m)
3. TDecompChol: Solve/TrnasSolve had incorrect
forward/
backward substitution function ! Tests
have been added to stressLinear to
check
Added NormlEqn solver which can handle
weights in data
4. TDecompSVD: Bug in Solve function ! index was
stop-
ping one too early . Was not noticed
by
stressLinear because the linear
equation
test wasnumerically to difficult
causing accuracy tolerances to be not
strict enough (was 0.005 , is now
DBL_EPSILON on a diagonal dominant
matrix)
Revision
8062 -
(
view)
(
download)
(
as text)
(
annotate)
-
[select for diffs]
Modified
Wed Jan 28 07:39:18 2004 UTC (10 years, 11 months ago) by
brun
Original Path:
trunk/matrix/inc/TDecompLU.h
File length: 2536 byte(s)
Diff to
previous 8013
From Eddy Offermann:
Reacting to Philippe's worry , I have removed
one TDecompLU constructor .
Answering George's Irwin question, I noticed
that I forgot to forward a change for the
inversion of symmetric matrices: I was
using the Cholesky decomp which wants the matrix
to be also positive definite
The patch rectifies these issues .
Revision
8013 -
(
view)
(
download)
(
as text)
(
annotate)
-
[select for diffs]
Added
Sun Jan 25 20:33:32 2004 UTC (11 years ago) by
brun
Original Path:
trunk/matrix/inc/TDecompLU.h
File length: 2588 byte(s)
New Linear Algebra package from Eddy Offermann.
This new package reimplements the previous classes TMatrix and TMatrixD.
The new classes should be back compatible with the previous version except
the function GetElements.
New classes have been introduced for symmetric matrices,
lazy matrices.
New algorithms (LU, SVD) have been introduced.
A new test suite test/stressLinear.cxx is introduced.
A complete description of this package will be posted in the coming days.
The classes are well documented in the implementation headers.
This form allows you to request diffs between any two revisions of this file.
For each of the two "sides" of the diff,
enter a numeric revision.
