Go to the documentation of this file.00001
00002
00003
00006
00007 #define WANT_MATH
00008
00009 #include "include.h"
00010
00011 #include "newmatap.h"
00012
00013 #ifdef use_namespace
00014 namespace NEWMAT {
00015 #endif
00016
00017
00018 #ifdef DO_REPORT
00019 #define REPORT { static ExeCounter ExeCount(__LINE__,21); ++ExeCount; }
00020 #else
00021 #define REPORT {}
00022 #endif
00023
00024 ReturnMatrix Helmert(int n, bool full)
00025 {
00026 REPORT
00027 Tracer et("Helmert ");
00028 if (n <= 0) Throw(ProgramException("Dimension <= 0 "));
00029 Matrix H;
00030
00031 if (full) H.resize(n,n); else H.resize(n-1,n);
00032 H = 0.0;
00033 for (int i = 1; i < n; ++i)
00034 {
00035 Real f = 1.0 / sqrt((Real)i * (i+1));
00036 H.submatrix(i,i,1,i) = -f; H(i,i+1) = f * i;
00037 }
00038 if (full) { H.row(n) = 1.0 / sqrt((Real)n); }
00039 H.release(); return H.for_return();
00040 }
00041
00042
00043
00044
00045
00046 ReturnMatrix Helmert(const ColumnVector& X, bool full)
00047 {
00048 REPORT
00049 Tracer et("Helmert * CV");
00050 int n = X.nrows();
00051 if (n == 0) Throw(ProgramException("X Vector of length 0", X));
00052 Real sum = 0.0; ColumnVector Y;
00053 if (full) Y.resize(n); else Y.resize(n-1);
00054 for (int i = 1; i < n; ++i)
00055 { sum += X(i); Y(i) = (i * X(i+1) - sum) / sqrt((Real)i * (i+1)); }
00056 if (full) { sum += X(n); Y(n) = sum / sqrt((Real)n); }
00057 Y.release(); return Y.for_return();
00058 }
00059
00060
00061 ReturnMatrix Helmert(int n, int j, bool full)
00062 {
00063 REPORT
00064 Tracer et("Helmert:single element ");
00065 if (n <= 0) Throw(ProgramException("X Vector of length <= 0"));
00066 if (j > n || j <= 0)
00067 Throw(ProgramException("Out of range element number "));
00068 ColumnVector Y; if (full) Y.resize(n); else Y.resize(n-1);
00069 Y = 0.0;
00070 if (j > 1) Y(j-1) = sqrt((Real)(j-1) / (Real)j);
00071 for (int i = j; i < n; ++i) Y(i) = - 1.0 / sqrt((Real)i * (i+1));
00072 if (full) Y(n) = 1.0 / sqrt((Real)n);
00073 Y.release(); return Y.for_return();
00074 }
00075
00076 ReturnMatrix Helmert_transpose(const ColumnVector& Y, bool full)
00077 {
00078 REPORT
00079 Tracer et("Helmert_transpose * CV ");
00080 int n = Y.nrows(); Real sum;
00081 if (full) sum = Y(n) / sqrt((Real)n); else { sum = 0.0; ++n; }
00082 ColumnVector X(n);
00083 for (int i = n-1; i > 0; --i)
00084 {
00085 Real Yi = Y(i) / sqrt((Real)i * (i+1));
00086 X(i+1) = i * Yi + sum; sum -= Yi;
00087 }
00088 X(1) = sum;
00089 X.release(); return X.for_return();
00090 }
00091
00092
00093 Real Helmert_transpose(const ColumnVector& Y, int j, bool full)
00094 {
00095 REPORT
00096 Tracer et("Helmert_transpose:single element ");
00097 int n = Y.nrows(); Real sum;
00098 if (full) sum = Y(n) / sqrt((Real)n); else { sum = 0.0; ++n; }
00099 if (j > n || j <= 0) Throw(IndexException(j, Y));
00100 for (int i = n-1; i > 0; --i)
00101 {
00102 Real Yi = Y(i) / sqrt((Real)i * (i+1));
00103 if (i+1 == j) return i * Yi + sum;
00104 sum -= Yi;
00105 }
00106 return sum;
00107 }
00108
00109 ReturnMatrix Helmert(const Matrix& X, bool full)
00110 {
00111 REPORT
00112 Tracer et("Helmert * Matrix");
00113 int m = X.nrows(); int n = X.ncols();
00114 if (m == 0) Throw(ProgramException("Matrix has 0 rows ", X));
00115 Matrix Y;
00116 if (full) Y.resize(m,n); else Y.resize(m-1, n);
00117 for (int j = 1; j <= n; ++j)
00118 {
00119 ColumnVector CV = X.Column(j);
00120 Y.Column(j) = Helmert(CV, full);
00121 }
00122 Y.release(); return Y.for_return();
00123 }
00124
00125 ReturnMatrix Helmert_transpose(const Matrix& Y, bool full)
00126 {
00127 REPORT
00128 Tracer et("Helmert_transpose * Matrix ");
00129 int m = Y.nrows(); int n = Y.ncols(); if (!full) ++m;
00130 Matrix X(m, n);
00131 for (int j = 1; j <= n; ++j)
00132 {
00133 ColumnVector CV = Y.Column(j);
00134 X.Column(j) = Helmert_transpose(CV, full);
00135 }
00136 X.release(); return X.for_return();
00137 }
00138
00139
00140
00141
00142 #ifdef use_namespace
00143 }
00144 #endif
00145
00146