00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00036 #ifndef MAT_UTILS_HEADER
00037 #define MAT_UTILS_HEADER
00038 #include "Interval.h"
00039 #include "matrix_proxy.h"
00040 namespace mat {
00041
00042 template<typename Tmatrix, typename Treal>
00043 struct DiffMatrix {
00044 typedef typename Tmatrix::VectorType VectorType;
00045 void getCols(SizesAndBlocks & colsCopy) const {
00046 A.getCols(colsCopy);
00047 }
00048 int get_nrows() const {
00049 assert( A.get_nrows() == B.get_nrows() );
00050 return A.get_nrows();
00051 }
00052 Treal frob() const {
00053 return Tmatrix::frob_diff(A, B);
00054 }
00055 void quickEuclBounds(Treal & euclLowerBound,
00056 Treal & euclUpperBound) const {
00057 Treal frobTmp = frob();
00058 euclLowerBound = frobTmp / template_blas_sqrt( (Treal)get_nrows() );
00059 euclUpperBound = frobTmp;
00060 }
00061
00062 Tmatrix const & A;
00063 Tmatrix const & B;
00064 DiffMatrix(Tmatrix const & A_, Tmatrix const & B_)
00065 : A(A_), B(B_) {}
00066 template<typename Tvector>
00067 void matVecProd(Tvector & y, Tvector const & x) const {
00068 Tvector tmp(y);
00069 tmp = (Treal)-1.0 * B * x;
00070 y = (Treal)1.0 * A * x;
00071 y += (Treal)1.0 * tmp;
00072 }
00073 };
00074
00075
00076
00077 template<typename Tmatrix, typename Treal>
00078 struct ATAMatrix {
00079 typedef typename Tmatrix::VectorType VectorType;
00080 Tmatrix const & A;
00081 explicit ATAMatrix(Tmatrix const & A_)
00082 : A(A_) {}
00083 void getCols(SizesAndBlocks & colsCopy) const {
00084 A.getRows(colsCopy);
00085 }
00086 void quickEuclBounds(Treal & euclLowerBound,
00087 Treal & euclUpperBound) const {
00088 Treal frobA = A.frob();
00089 euclLowerBound = 0;
00090 euclUpperBound = frobA * frobA;
00091 }
00092
00093
00094 template<typename Tvector>
00095 void matVecProd(Tvector & y, Tvector const & x) const {
00096 y = x;
00097 y = A * y;
00098 y = transpose(A) * y;
00099 }
00100
00101 int get_nrows() const { return A.get_ncols(); }
00102 };
00103
00104
00105 template<typename Tmatrix, typename Tmatrix2, typename Treal>
00106 struct TripleMatrix {
00107 typedef typename Tmatrix::VectorType VectorType;
00108 void getCols(SizesAndBlocks & colsCopy) const {
00109 A.getCols(colsCopy);
00110 }
00111 int get_nrows() const {
00112 assert( A.get_nrows() == Z.get_nrows() );
00113 return A.get_nrows();
00114 }
00115 void quickEuclBounds(Treal & euclLowerBound,
00116 Treal & euclUpperBound) const {
00117 Treal frobA = A.frob();
00118 Treal frobZ = Z.frob();
00119 euclLowerBound = 0;
00120 euclUpperBound = frobA * frobZ * frobZ;
00121 }
00122
00123 Tmatrix const & A;
00124 Tmatrix2 const & Z;
00125 TripleMatrix(Tmatrix const & A_, Tmatrix2 const & Z_)
00126 : A(A_), Z(Z_) {}
00127 void matVecProd(VectorType & y, VectorType const & x) const {
00128 VectorType tmp(x);
00129 tmp = Z * tmp;
00130 y = (Treal)1.0 * A * tmp;
00131 y = transpose(Z) * y;
00132 }
00133 };
00134
00135
00136 template<typename Tmatrix, typename Tmatrix2, typename Treal>
00137 struct CongrTransErrorMatrix {
00138 typedef typename Tmatrix::VectorType VectorType;
00139 void getCols(SizesAndBlocks & colsCopy) const {
00140 E.getRows(colsCopy);
00141 }
00142 int get_nrows() const {
00143 return E.get_ncols();
00144 }
00145 void quickEuclBounds(Treal & euclLowerBound,
00146 Treal & euclUpperBound) const {
00147 Treal frobA = A.frob();
00148 Treal frobZ = Zt.frob();
00149 Treal frobE = E.frob();
00150 euclLowerBound = 0;
00151 euclUpperBound = frobA * frobE * frobE + 2 * frobA * frobE * frobZ;
00152 }
00153
00154 Tmatrix const & A;
00155 Tmatrix2 const & Zt;
00156 Tmatrix2 const & E;
00157
00158 CongrTransErrorMatrix(Tmatrix const & A_,
00159 Tmatrix2 const & Z_,
00160 Tmatrix2 const & E_)
00161 : A(A_), Zt(Z_), E(E_) {}
00162 void matVecProd(VectorType & y, VectorType const & x) const {
00163
00164 VectorType tmp(x);
00165 tmp = E * tmp;
00166 y = (Treal)-1.0 * A * tmp;
00167 y = transpose(E) * y;
00168
00169 VectorType tmp1;
00170 tmp = x;
00171 tmp = Zt * tmp;
00172 tmp1 = (Treal)1.0 * A * tmp;
00173 tmp1 = transpose(E) * tmp1;
00174 y += (Treal)1.0 * tmp1;
00175
00176 tmp = x;
00177 tmp = E * tmp;
00178 tmp1 = (Treal)1.0 * A * tmp;
00179 tmp1 = transpose(Zt) * tmp1;
00180 y += (Treal)1.0 * tmp1;
00181 }
00182 };
00183
00184
00185
00186 }
00187 #endif