jama_qr.h Source File
Main Page Namespace List Compound List File List Compound Members
Copyright By PowCoder代写 加微信 powcoder
Go to the documentation of this file.00001 #ifndef JAMA_QR_H
00002 #define JAMA_QR_H
00004 #include “tnt_array1d.h”
00005 #include “tnt_array2d.h”
00006 #include “tnt_math_utils.h”
00008 namespace JAMA
00034 template
00035 class QR {
00042 TNT::Array2D
00048 int m, n;
00053 TNT::Array1D
00056 public:
00063 QR(const TNT::Array2D
00064 {
00065 QR_ = A.copy();
00066 m = A.dim1();
00067 n = A.dim2();
00068 Rdiag = TNT::Array1D
00069 int i=0, j=0, k=0;
00071 // Main loop.
00072 for (k = 0; k < n; k++) {
00073 // Compute 2-norm of k-th column without under/overflow.
00074 Real nrm = 0;
00075 for (i = k; i < m; i++) {
00076 nrm = hypot(nrm,QR_[i][k]);
00077 }
00079 if (nrm != 0.0) {
00080 // Form k-th Householder vector.
00081 if (QR_[k][k] < 0) {
00082 nrm = -nrm;
00083 }
00084 for (i = k; i < m; i++) {
00085 QR_[i][k] /= nrm;
00086 }
00087 QR_[k][k] += 1.0;
00089 // Apply transformation to remaining columns.
00090 for (j = k+1; j < n; j++) {
00091 Real s = 0.0;
00092 for (i = k; i < m; i++) {
00093 s += QR_[i][k]*QR_[i][j];
00094 }
00095 s = -s/QR_[k][k];
00096 for (i = k; i < m; i++) {
00097 QR_[i][j] += s*QR_[i][k];
00098 }
00099 }
00100 }
00101 Rdiag[k] = -nrm;
00102 }
00103 }
00111 int isFullRank() const
00112 {
00113 for (int j = 0; j < n; j++)
00114 {
00115 if (Rdiag[j] == 0)
00116 return 0;
00117 }
00118 return 1;
00119 }
00130 TNT::Array2D
00131 {
00132 TNT::Array2D
00134 /* note: H is completely filled in by algorithm, so
00135 initializaiton of H is not necessary.
00136 */
00137 for (int i = 0; i < m; i++)
00138 {
00139 for (int j = 0; j < n; j++)
00140 {
00141 if (i >= j) {
00142 H[i][j] = QR_[i][j];
00143 } else {
00144 H[i][j] = 0.0;
00145 }
00146 }
00147 }
00148 return H;
00149 }
00157 TNT::Array2D
00158 {
00159 TNT::Array2D
00160 for (int i = 0; i < n; i++) {
00161 for (int j = 0; j < n; j++) {
00162 if (i < j) {
00163 R[i][j] = QR_[i][j];
00164 } else if (i == j) {
00165 R[i][j] = Rdiag[i];
00166 } else {
00167 R[i][j] = 0.0;
00168 }
00169 }
00170 }
00171 return R;
00172 }
00183 TNT::Array2D
00184 {
00185 int i=0, j=0, k=0;
00187 TNT::Array2D
00188 for (k = n-1; k >= 0; k–) {
00189 for (i = 0; i < m; i++) {
00190 Q[i][k] = 0.0;
00191 }
00192 Q[k][k] = 1.0;
00193 for (j = k; j < n; j++) {
00194 if (QR_[k][k] != 0) {
00195 Real s = 0.0;
00196 for (i = k; i < m; i++) {
00197 s += QR_[i][k]*Q[i][j];
00198 }
00199 s = -s/QR_[k][k];
00200 for (i = k; i < m; i++) {
00201 Q[i][j] += s*QR_[i][k];
00202 }
00203 }
00204 }
00205 }
00206 return Q;
00207 }
00217 TNT::Array1D
00218 {
00219 if (b.dim1() != m) /* arrays must be conformant */
00220 return TNT::Array1D
00222 if ( !isFullRank() ) /* matrix is rank deficient */
00223 {
00224 return TNT::Array1D
00225 }
00227 TNT::Array1D
00228 int i=0, j=0, k=0;
00230 // Compute Y = transpose(Q)*b
00231 for (k = 0; k < n; k++)
00232 {
00233 Real s = 0.0;
00234 for (i = k; i < m; i++)
00235 {
00236 s += QR_[i][k]*x[i];
00237 }
00238 s = -s/QR_[k][k];
00239 for (i = k; i < m; i++)
00240 {
00241 x[i] += s*QR_[i][k];
00242 }
00243 }
00244 // Solve R*X = Y;
00245 for (k = n-1; k >= 0; k–)
00246 {
00247 x[k] /= Rdiag[k];
00248 for (i = 0; i < k; i++) {
00249 x[i] -= x[k]*QR_[i][k];
00250 }
00251 }
00254 /* return n x nx portion of X */
00255 TNT::Array1D
00256 for (i=0; i
00270 {
00271 if (B.dim1() != m) /* arrays must be conformant */
00272 return TNT::Array2D
00274 if ( !isFullRank() ) /* matrix is rank deficient */
00275 {
00276 return TNT::Array2D
00277 }
00279 int nx = B.dim2();
00280 TNT::Array2D
00281 int i=0, j=0, k=0;
00283 // Compute Y = transpose(Q)*B
00284 for (k = 0; k < n; k++) {
00285 for (j = 0; j < nx; j++) {
00286 Real s = 0.0;
00287 for (i = k; i < m; i++) {
00288 s += QR_[i][k]*X[i][j];
00289 }
00290 s = -s/QR_[k][k];
00291 for (i = k; i < m; i++) {
00292 X[i][j] += s*QR_[i][k];
00293 }
00294 }
00295 }
00296 // Solve R*X = Y;
00297 for (k = n-1; k >= 0; k–) {
00298 for (j = 0; j < nx; j++) {
00299 X[k][j] /= Rdiag[k];
00300 }
00301 for (i = 0; i < k; i++) {
00302 for (j = 0; j < nx; j++) {
00303 X[i][j] -= X[k][j]*QR_[i][k];
00304 }
00305 }
00306 }
00309 /* return n x nx portion of X */
00310 TNT::Array2D
00311 for (i=0; i