685{
686 if(ncol != nrow)
687 {
688 error(
"G4ErrorMatrix::invert: G4ErrorMatrix is not NxN");
689 }
690
693 if(!ir)
694 ir =
new G4int[max_array + 1];
695
696 if(ncol > max_array)
697 {
698 delete[] ir;
699 max_array = nrow;
700 ir =
new G4int[max_array + 1];
701 }
705 switch(nrow)
706 {
707 case 3:
708 G4double c11, c12, c13, c21, c22, c23, c31, c32, c33;
709 ifail = 0;
710 c11 = (*(m.begin() + 4)) * (*(m.begin() + 8)) -
711 (*(m.begin() + 5)) * (*(m.begin() + 7));
712 c12 = (*(m.begin() + 5)) * (*(m.begin() + 6)) -
713 (*(m.begin() + 3)) * (*(m.begin() + 8));
714 c13 = (*(m.begin() + 3)) * (*(m.begin() + 7)) -
715 (*(m.begin() + 4)) * (*(m.begin() + 6));
716 c21 = (*(m.begin() + 7)) * (*(m.begin() + 2)) -
717 (*(m.begin() + 8)) * (*(m.begin() + 1));
718 c22 = (*(m.begin() + 8)) * (*m.begin()) -
719 (*(m.begin() + 6)) * (*(m.begin() + 2));
720 c23 = (*(m.begin() + 6)) * (*(m.begin() + 1)) -
721 (*(m.begin() + 7)) * (*m.begin());
722 c31 = (*(m.begin() + 1)) * (*(m.begin() + 5)) -
723 (*(m.begin() + 2)) * (*(m.begin() + 4));
724 c32 = (*(m.begin() + 2)) * (*(m.begin() + 3)) -
725 (*m.begin()) * (*(m.begin() + 5));
726 c33 = (*m.begin()) * (*(m.begin() + 4)) -
727 (*(m.begin() + 1)) * (*(m.begin() + 3));
728 t1 = std::fabs(*m.begin());
729 t2 = std::fabs(*(m.begin() + 3));
730 t3 = std::fabs(*(m.begin() + 6));
731 if(t1 >= t2)
732 {
733 if(t3 >= t1)
734 {
735 temp = *(m.begin() + 6);
736 det = c23 * c12 - c22 * c13;
737 }
738 else
739 {
740 temp = *(m.begin());
741 det = c22 * c33 - c23 * c32;
742 }
743 }
744 else if(t3 >= t2)
745 {
746 temp = *(m.begin() + 6);
747 det = c23 * c12 - c22 * c13;
748 }
749 else
750 {
751 temp = *(m.begin() + 3);
752 det = c13 * c32 - c12 * c33;
753 }
754 if(det == 0)
755 {
756 ierr = 1;
757 return;
758 }
759 {
760 ss = temp / det;
762 *(mq++) = ss * c11;
763 *(mq++) = ss * c21;
764 *(mq++) = ss * c31;
765 *(mq++) = ss * c12;
766 *(mq++) = ss * c22;
767 *(mq++) = ss * c32;
768 *(mq++) = ss * c13;
769 *(mq++) = ss * c23;
770 *(mq) = ss * c33;
771 }
772 break;
773 case 2:
774 ifail = 0;
775 det = (*m.begin()) * (*(m.begin() + 3)) -
776 (*(m.begin() + 1)) * (*(m.begin() + 2));
777 if(det == 0)
778 {
779 ierr = 1;
780 return;
781 }
782 ss = 1.0 / det;
783 temp = ss * (*(m.begin() + 3));
784 *(m.begin() + 1) *= -ss;
785 *(m.begin() + 2) *= -ss;
786 *(m.begin() + 3) = ss * (*m.begin());
787 *(m.begin()) = temp;
788 break;
789 case 1:
790 ifail = 0;
791 if((*(m.begin())) == 0)
792 {
793 ierr = 1;
794 return;
795 }
796 *(m.begin()) = 1.0 / (*(m.begin()));
797 break;
798 case 4:
800 return;
801 case 5:
803 return;
804 case 6:
806 return;
807 default:
808 ifail = dfact_matrix(det, ir);
809 if(ifail)
810 {
811 ierr = 1;
812 return;
813 }
814 dfinv_matrix(ir);
815 break;
816 }
817 ierr = 0;
818 return;
819}
virtual void invertHaywood4(G4int &ierr)
virtual void invertHaywood5(G4int &ierr)
virtual void invertHaywood6(G4int &ierr)