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