LinAlg: detsymm.cpp Source File

00001 //===========================================================================
00043 //===========================================================================
00044 
00045 #include <cmath>
00046 #include <SharkDefs.h>
00047 #include <LinAlg/LinAlg.h>
00048 
00049 //===========================================================================
00093 double detsymm
00094 (
00095     Array2D< double >& amatA,
00096     Array2D< double >& vmatA,
00097     Array  < double >& dvecA
00098 )
00099 {
00100     SIZE_CHECK(amatA.dim(0) == amatA.dim(1))
00101 
00102     unsigned n = amatA.dim(0);
00103     unsigned i;
00104     unsigned j;
00105     double   det;
00106 
00107     // Fill upper triangular matrix:
00108     for (i = 0; i + 1 < n; i++) {
00109         for (j = i + 1; j < n; j++) {
00110             amatA(i, j) = amatA(j, i);
00111         }
00112     }
00113 
00114     // Calculate eigenvalues and use intermediate results stored in hmatA
00115     //eigensymm( amatA, vmatA, dvecA );
00116     Array< double > hmatA(amatA);
00117     eigensymm_intermediate(amatA, hmatA, vmatA, dvecA);
00118 
00119     for (i = 0; i + 1 < n; i++) {
00120         for (j = i + 1; j < n; j++) {
00121             hmatA(j, i) = hmatA(i, j);
00122             amatA(j, i) = hmatA(j, i);
00123             amatA(i, j) = hmatA(i, j);
00124         }
00125     }
00126 
00127     // Calculate determinate as product of eigenvalues:
00128     for (i = 0, det = 1; i < n; det *= dvecA(i++));
00129 
00130     return det;
00131 }
00132 
00133 
00134 //===========================================================================
00175 double logdetsymm
00176 (
00177     Array2D< double >& amatA,
00178     Array2D< double >& vmatA,
00179     Array  < double >& dvecA
00180 )
00181 {
00182     SIZE_CHECK(amatA.dim(0) == amatA.dim(1))
00183 
00184     unsigned n = amatA.dim(0);
00185     unsigned i;
00186     unsigned j;
00187     double   det = 0.;
00188 
00189     // Fill upper triangular matrix:
00190     for (i = 0; i + 1 < n; i++) {
00191         for (j = i + 1; j < n; j++) {
00192             amatA(i, j) = amatA(j, i);
00193         }
00194     }
00195 
00196     // Calculate eigenvalues and use intermediate results stored in hmatA
00197     //eigensymm( amatA, vmatA, dvecA );
00198     Array< double > hmatA(amatA);
00199     eigensymm_intermediate(amatA, hmatA, vmatA, dvecA);
00200 
00201     for (i = 0; i + 1 < n; i++) {
00202         for (j = i + 1; j < n; j++) {
00203             hmatA(j, i) = hmatA(i, j);
00204             amatA(j, i) = hmatA(j, i);
00205             amatA(i, j) = hmatA(i, j);
00206         }
00207     }
00208 
00209     // Calculate determinate as product of eigenvalues:
00210     for (i = 0; i < n; det += log(dvecA(i++)));
00211 
00212     return det;
00213 }