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 #ifndef QR_H
00026 #define QR_H
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050 #include <cmath>
00051 #include "tntmath.h"
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081 namespace SPUC
00082 {
00083
00084 template <class MaTRiX, class Vector>
00085 int QR_factor(MaTRiX &A, Vector& C, Vector &D)
00086 {
00087 assert(A.lbound() == 1);
00088 assert(C.lbound() == 1);
00089 assert(D.lbound() == 1);
00090
00091 int M = A.dim1();
00092 int N = A.dim2();
00093
00094 assert(M == N);
00095
00096 int i,j,k;
00097 typename MaTRiX::element_type eta, sigma, sum;
00098
00099
00100
00101 if (N != C.size()) C.newsize(N);
00102 if (N != D.size()) D.newsize(N);
00103
00104 for (k=1; k<N; k++)
00105 {
00106
00107
00108 eta = 0;
00109 for (i=k; i<=N; i++)
00110 {
00111 double absA = fabs(A(i,k));
00112 eta = ( absA > eta ? absA : eta );
00113 }
00114
00115 if (eta == 0)
00116 {
00117 cerr << "QR: k=" << k << "\n";
00118 return 1;
00119 }
00120
00121
00122
00123 for(i=k; i<=N; i++)
00124 A(i,k) = A(i,k) / eta;
00125
00126 sum = 0;
00127 for (i=k; i<=N; i++)
00128 sum = sum + A(i,k)*A(i,k);
00129 sigma = sign(A(k,k)) * root(sum);
00130
00131
00132 A(k,k) = A(k,k) + sigma;
00133 C(k) = sigma * A(k,k);
00134 D(k) = -eta * sigma;
00135
00136 for (j=k+1; j<=N; j++)
00137 {
00138 sum = 0;
00139 for (i=k; i<=N; i++)
00140 sum = sum + A(i,k)*A(i,j);
00141 sum = sum / C(k);
00142
00143 for (i=k; i<=N; i++)
00144 A(i,j) = A(i,j) - sum * A(i,k);
00145 }
00146
00147 D(N) = A(N,N);
00148 }
00149
00150 return 0;
00151 }
00152
00153
00154
00155
00156 template <class MaTRiX, class Vector>
00157 int R_solve(const MaTRiX &A, Vector &D, Vector &b)
00158 {
00159 assert(A.lbound() == 1);
00160 assert(D.lbound() == 1);
00161 assert(b.lbound() == 1);
00162
00163 int i,j;
00164 int N = A.dim1();
00165
00166 assert(N == A.dim2());
00167 assert(N == D.dim());
00168 assert(N == b.dim());
00169
00170 typename MaTRiX::element_type sum;
00171
00172 if (D(N) == 0)
00173 return 1;
00174
00175 b(N) = b(N) /
00176 D(N);
00177
00178 for (i=N-1; i>=1; i--)
00179 {
00180 if (D(i) == 0)
00181 return 1;
00182 sum = 0;
00183 for (j=i+1; j<=N; j++)
00184 sum = sum + A(i,j)*b(j);
00185 b(i) = ( b(i) - sum ) /
00186 D(i);
00187 }
00188
00189 return 0;
00190 }
00191
00192
00193 template <class MaTRiX, class Vector>
00194 int QR_solve(const MaTRiX &A, const Vector &c, Vector &d,
00195 Vector &b)
00196 {
00197 assert(A.lbound() == 1);
00198 assert(c.lbound() == 1);
00199 assert(d.lbound() == 1);
00200
00201 int N=A.dim1();
00202
00203 assert(N == A.dim2());
00204 assert(N == c.dim());
00205 assert(N == d.dim());
00206 assert(N == b.dim());
00207
00208 int i,j;
00209 typename MaTRiX::element_type sum, tau;
00210
00211 for (j=1; j<N; j++)
00212 {
00213
00214 sum = 0;
00215 for (i=j; i<=N; i++)
00216 sum = sum + A(i,j)*b(i);
00217 if (c(j) == 0)
00218 return 1;
00219 tau = sum / c(j);
00220 for (i=j; i<=N; i++)
00221 b(i) = b(i) - tau * A(i,j);
00222 }
00223 return R_solve(A, d, b);
00224 }
00225
00226 }
00227
00228 #endif
00229