Changeset 8901 in ntrip for trunk/BNC/newmat/hholder.cpp
- Timestamp:
- Mar 18, 2020, 11:06:13 AM (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/BNC/newmat/hholder.cpp
r2013 r8901 335 335 } 336 336 337 // Following previous transformation, 338 // now apply the same orthogonal transformation to (MX & MU) 339 // Need the X Matrix but not the U. 340 // Not optimised for accessing consecutive memory 341 342 void updateQRZ(const Matrix& X, Matrix& MX, Matrix& MU) 343 { 344 REPORT 345 Tracer et("updateQRZ(2)"); 346 int s = X.Ncols(); int n = X.Nrows(); 347 if (n != MX.Nrows()) 348 Throw(ProgramException("Incompatible dimensions",X,MX)); 349 if (s != MU.Nrows()) 350 Throw(ProgramException("Incompatible dimensions",X,MU)); 351 int t = MX.Ncols(); 352 if (t != MU.Ncols()) 353 Throw(ProgramException("Incompatible dimensions",MX,MU)); 354 355 if (s == 0) return; 356 357 const Real* xi0 = X.data(); Real* mx = MX.data(); Real* muj = MU.data(); 358 for (int i=1; i<=s; ++i) 359 { 360 Real sum = 0.0; 361 { 362 const Real* xi=xi0; int k=n; 363 while(k--) { sum += square(*xi); xi+= s;} 364 } 365 Real a0 = sqrt(2.0 - sum); Real* mxj0 = mx; 366 for (int j=1; j<=t; ++j) 367 { 368 Real sum = 0.0; 369 const Real* xi=xi0; Real* mxj=mxj0; int k=n; 370 while(--k) { sum += *xi * *mxj; xi += s; mxj += t; } 371 sum += *xi * *mxj; // last line of loop 372 sum += a0 * *muj; 373 xi=xi0; mxj=mxj0; k=n; 374 while(--k) { *mxj -= sum * *xi; xi += s; mxj += t; } 375 *mxj -= sum * *xi; // last line of loop 376 *muj -= sum * a0; ++mxj0; ++muj; 377 } 378 ++xi0; 379 } 380 } 381 382 383 384 // same thing as updateQRZ(Matrix& X, UpperTriangularMatrix& U) 385 // except that X is upper triangular 386 // contents of X are destroyed - results are in U 387 // assume we can access efficiently by columns 388 // e.g. X and U will fit in cache memory 389 390 void updateQRZ(UpperTriangularMatrix& X, UpperTriangularMatrix& U) 391 { 392 REPORT 393 Tracer et("updateQRZ(3)"); 394 int s = X.Ncols(); 395 if (s != U.Ncols()) 396 Throw(ProgramException("Incompatible dimensions",X,U)); 397 if (s == 0) return; 398 Real* xi0 = X.data(); Real* u = U.data(); 399 for (int i=1; i<=s; ++i) 400 { 401 Real r = *u; Real sum = 0.0; 402 { 403 Real* xi=xi0; int k=i; int l=s; 404 while(k--) { sum += square(*xi); xi+= --l;} 405 } 406 sum = sqrt(sum + square(r)); 407 if (sum == 0.0) { REPORT X.column(i) = 0.0; *u = 0.0; } 408 else 409 { 410 Real frs = fabs(r) + sum; 411 Real a0 = sqrt(frs / sum); Real alpha = a0 / frs; 412 if (r <= 0) { REPORT *u = sum; alpha = -alpha; } 413 else { REPORT *u = -sum; } 414 { 415 Real* xj0=xi0; int k=i; int l=s; 416 while(k--) { *xj0 *= alpha; --l; xj0 += l;} 417 } 418 Real* xj0=xi0; Real* uj=u; 419 for (int j=i+1; j<=s; ++j) 420 { 421 Real sum = 0.0; ++xj0; ++uj; 422 Real* xi=xi0; Real* xj=xj0; int k=i; int l=s; 423 while(k--) { sum += *xi * *xj; --l; xi += l; xj += l; } 424 sum += a0 * *uj; 425 xi=xi0; xj=xj0; k=i; l=s; 426 while(k--) { *xj -= sum * *xi; --l; xi += l; xj += l; } 427 *uj -= sum * a0; 428 } 429 } 430 ++xi0; u += s-i+1; 431 } 432 } 433 434 // Following previous transformation, 435 // now apply the same orthogonal transformation to (MX & MU) 436 // Need the X UpperTriangularMatrix but not the U. 437 438 void updateQRZ(const UpperTriangularMatrix& X, Matrix& MX, Matrix& MU) 439 { 440 REPORT 441 Tracer et("updateQRZ(4)"); 442 int s = X.Ncols(); 443 if (s != MX.Nrows()) 444 Throw(ProgramException("Incompatible dimensions",X,MX)); 445 if (s != MU.Nrows()) 446 Throw(ProgramException("Incompatible dimensions",X,MU)); 447 int t = MX.Ncols(); 448 if (t != MU.Ncols()) 449 Throw(ProgramException("Incompatible dimensions",MX,MU)); 450 if (s == 0) return; 451 452 const Real* xi0 = X.data(); Real* mx = MX.data(); Real* muj = MU.data(); 453 for (int i=1; i<=s; ++i) 454 { 455 Real sum = 0.0; 456 { 457 const Real* xi=xi0; int k=i; int l=s; 458 while(k--) { sum += square(*xi); xi+= --l;} 459 } 460 Real a0 = sqrt(2.0 - sum); Real* mxj0 = mx; 461 for (int j=1; j<=t; ++j) 462 { 463 Real sum = 0.0; 464 const Real* xi=xi0; Real* mxj=mxj0; int k=i; int l=s; 465 while(--k) { sum += *xi * *mxj; --l; xi += l; mxj += t; } 466 sum += *xi * *mxj; // last line of loop 467 sum += a0 * *muj; 468 xi=xi0; mxj=mxj0; k=i; l=s; 469 while(--k) { *mxj -= sum * *xi; --l; xi += l; mxj += t; } 470 *mxj -= sum * *xi; // last line of loop 471 *muj -= sum * a0; ++mxj0; ++muj; 472 } 473 ++xi0; 474 } 475 } 476 477 478 479 480 337 481 // Matrix A's first n columns are orthonormal 338 482 // so A.Columns(1,n).t() * A.Columns(1,n) is the identity matrix.
Note:
See TracChangeset
for help on using the changeset viewer.