368 {
369
370
371
372 cout << std::setiosflags(std::ios::fixed) << std::setw(10);
377
378
379
380
382
383
384
390
391
392
396#if defined(PRINTOUT)
397 cout << a << endl;
398#endif
401#if defined(PRINTOUT)
402 cout << b << endl;
403#endif
404
411#if defined(PRINTOUT)
412 cout << "nraw=" << b.num_row() << " ncol=" << b.num_col() << endl;
413#endif
414#if defined(PRINTOUT)
415 cout << "b(1,1)=" << b(1,1) << " b(2,1)=" << b(2,1) << endl;
416#endif
417 b(2,3) = 1.0;
418
419 b /= 3.0;
420 b *= 6.0;
421
425
426 d += e;
427 d += sm;
428 d += dm;
429 ddmv += vvm;
430
431 a -= e;
432 a -= sm;
433 a -= dm;
434 ddmv -= vvm;
435
436 d = sm;
437 d = dm;
438 d = a;
439 ddmv = vvm;
440
442
443 a = d.T();
444
445 e = b.sub(1,2,3,2);
446
447 b.sub(1,3,e);
448
449 swap(a,d);
450
451 m = d * a;
452 m = d * sm;
453 m = d * dm;
454 m = sm * d;
455#if defined(PRINTOUT)
456 cout << "Dm=" << dm << " d = " << d << endl;
457#endif
458 m = dm * d;
459 m = sm * sm;
460 m = dm * dm;
461
462
463
464
472
474#if defined(PRINTOUT)
475 cout << "nraw=" << sr.num_row() << " ncol=" << sr.num_col() << endl;
476#endif
477#if defined(PRINTOUT)
478 cout << "sr(1,1)=" << sr(1,1) << " sr(2,1)=" << sr(2,1) << endl;
479#endif
480
481 sr(4,3) = r();
482 sr.fast(3,1) = r();
483
485#if defined(PRINTOUT)
486 cout << "d=" << d << "s=" << s << endl;
487#endif
489
490 ss *= 0.5;
491 ss /= 0.2;
492
494 ss += si;
495 ss += ddds;
496 ss -= sr;
497 ss -= ddds;
498 ss = ddds;
499
500 s = sr;
505
508
510
512
515
516 m = s * sr;
517 m = s * m;
518 m = m * s;
519
522
523
524
525
526
533
535#if defined(PRINTOUT)
536 cout << "Nr=" << di.num_row() << " Nc=" << dz.num_col() << endl;
537#endif
538#if defined(PRINTOUT)
539 cout << "dr(1,1)=" << dr(1,1) << endl;
540#endif
541
542 dr(4,4) = r();
543 dr.fast(2,2) = r();
544
548
549 dr *= 3.0;
550 dr /= 3.0;
551
552 dr += dt;
553 dr -= di;
554
555 dt = dz;
556
559
562#if defined(PRINTOUT)
564#endif
565
568
569
570
571
572
579#if defined(PRINTOUT)
580 cout << "vz=" << vz << "vi=" << vi << endl;
581#endif
582
586#if defined(PRINTOUT)
587 cout << "Nr=" << vr.num_row() << endl;
588#endif
589#if defined(PRINTOUT)
590 cout << "vr(3)=" << vr(3) << endl;
591#endif
592
593 vr(2) = r();
594
595 vi *= 2.5;
596 vi /= 2.5;
597
601 vi += vr;
602
606 vi -= vc;
607
611 vt = vc;
612
616#if defined(PRINTOUT)
617 cout <<
"Normsq=" << vc.
normsq() <<
"Norm=" << vc.norm() << endl;
618#endif
619
620 m = vc.T();
621
622 swap(vc,vr);
623
624 vt = sr * vr;
625 vt = dr * vr;
626 vt = mt * vr;
628
629
630
631
632
634 s = 3.5 * ss;
635 dt = 3.5 * dr;
636 vt = 3.5 * vr;
637
639 s = ss * 3.6;
640 dt = dr * 3.6;
641 vt = vr * 3.6;
642
644 s = ss / 3.6;
645 dt = dr / 3.6;
646 vt = vr / 3.6;
647
650 m = d + e;
652 m = d + sr;
653 m = d + dr;
654
655 m = sr + d;
656 m = dr + d;
657
658 s = sr + si;
659 dt = dr + di;
660
661 s = dr + sr;
662 s = sr + dr;
663
665
666 v = e + vr;
667
668 v = vr + e;
669
670 v = vr + vz;
671
672
675 m = d - e;
678 m = d - sr;
679 m = d - dr;
680
681 m = sr - d;
682 m = dr - d;
683
684 s = sr - si;
685 dt = dr - di;
686
687 s = dr - sr;
688 s = sr - dr;
689
691
692 v = e - vr;
693
694 v = vr - e;
695
696 v = vr - vz;
697
698
699
700
701
703
704 int inv;
705 a = m.inverse(inv);
706#if defined(PRINTOUT)
707 cout << "Inv=" << inv << endl;
708#endif
709 if(inv==0) {
710#if defined(PRINTOUT)
712 cout <<
"a*m=" << am.apply(
absf) <<
"m*a=" << ma.apply(
absf) << endl;
713#endif
714 }
715
717
719#if defined(PRINTOUT)
720 cout << "m*vt=" << m*vt << "v=" << v << endl;
721#endif
722
723
725
726 s = ss.inverse(inv);
727#if defined(PRINTOUT)
728 cout << "Inv=" << inv << endl;
729#endif
730 if(inv==0) {
731#if defined(PRINTOUT)
733 cout << sss.apply(
absf) << ssss.apply(
absf) << endl;
734#endif
735 }
736
737#if defined(PRINTOUT)
738 cout << "Before diag:ss=" << ss << endl;
739 s = ss;
740#endif
742#if defined(PRINTOUT)
743 cout << "m=" << m << endl;
744#endif
745#if defined(PRINTOUT)
746 cout << "ss=" << ss << endl;
747 cout << "Sim" << ss.similarity(m) << endl;
749 cout <<
"Diff" << diff.apply(
absf) << endl;
750#endif
751
754#if defined(PRINTOUT)
756 cout << am.apply(
absf) << ma.apply(
absf) << endl;
757#endif
758#if defined(PRINTOUT)
759 cout <<
"Norm 1 =" <<
norm1(m) << endl;
760#endif
761#if defined(PRINTOUT)
762 cout <<
"Norm 2 =" <<
norm(m) << endl;
763#endif
764#if defined(PRINTOUT)
766#endif
767
769#if defined(PRINTOUT)
770 cout << "m=" << m << endl;
771#endif
773#if defined(PRINTOUT)
774 cout << "a=" << a << "m=" << m << endl;
775#endif
776#if defined(PRINTOUT)
777 cout << "a*m=" << a*m << endl;
778#endif
779
783#if defined(PRINTOUT)
784 cout << "v=" << v << " m*vt=" << m*vt << endl;
785#endif
786
790#if defined(PRINTOUT)
791 cout << "v=" << a << " m*b=" << m*b << endl;
792#endif
793#if defined(PRINTOUT)
794 cout << "Test was successful" << endl;
795#endif
796
797
798
799
800
801
805 int i;
806 for(i=0;i<10;i++) {
808 p[i] *= 2.0;
809 p[i](1) -= 1;
810 p[i](2) -= 1;
811 p[i](3) -= 1;
812 }
814 for(i=0;i<10;i++) {
815 double psq = p[i].
normsq();
816 sp(1,1) += psq - p[i](1)*p[i](1);
817 sp(2,2) += psq - p[i](2)*p[i](2);
818 sp(3,3) += psq - p[i](3)*p[i](3);
819 sp(2,1) -= p[i](2)*p[i](1);
820 sp(3,1) -= p[i](3)*p[i](1);
821 sp(3,2) -= p[i](3)*p[i](2);
822 }
824 cout << "sp=" << sp << " spd=" << spd << endl;
827 cout <<
"condition = " <<
condition(spps) << endl;
830 cout <<
"spps=" << spps <<
"sppss - sim = " << sppssdiff.apply(
absf)
831 << endl;
832 return 0;
833}
void assign(const HepMatrix &hm2)
HepSymMatrix similarityT(const HepMatrix &hm1) const
HepDiagMatrix apply(double(*f)(double, int, int)) const
void assign(const HepMatrix &hm2)
HepSymMatrix apply(double(*f)(double, int, int)) const
HepSymMatrix similarityT(const HepMatrix &hm1) const
HepSymMatrix sub(int min_row, int max_row) const
HepSymMatrix similarity(const HepMatrix &hm1) const
HepVector sub(int min_row, int max_row) const
HepVector apply(double(*f)(double, int)) const
HepVector qr_solve(const HepMatrix &A, const HepVector &b)
HepSymMatrix vT_times_v(const HepVector &v)
double norm(const HepGenMatrix &m)
HepMatrix qr_inverse(const HepMatrix &A)
void qr_decomp(HepMatrix *A, HepMatrix *hsm)
double norm1(const HepGenMatrix &m)
double norm_infinity(const HepGenMatrix &m)
double condition(const HepSymMatrix &m)
HepMatrix diagonalize(HepSymMatrix *s)
HepVector solve(const HepMatrix &, const HepVector &)
double absf(double f, int, int)
double neg(double f, int, int)
double negv(double f, int)
int matrix_test1(const HepGenMatrix &m)