Garfield++ 4.0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
Functions
Garfield::Numerics::CERNLIB Namespace Reference

Linear algebra routines from CERNLIB. More...

Functions

int deqn (const int n, std::vector< std::vector< double > > &a, std::vector< double > &b)
 
int deqinv (const int n, std::vector< std::vector< double > > &a, std::vector< double > &b)
 Replaces b by the solution x of Ax = b, and replace A by its inverse.
 
void dfact (const int n, std::vector< std::vector< double > > &a, std::vector< int > &ir, int &ifail, double &det, int &jfail)
 
void dfeqn (const int n, std::vector< std::vector< double > > &a, std::vector< int > &ir, std::vector< double > &b)
 
void dfinv (const int n, std::vector< std::vector< double > > &a, std::vector< int > &ir)
 
int dinv (const int n, std::vector< std::vector< double > > &a)
 Replace square matrix A by its inverse.
 
void cfact (const int n, std::vector< std::vector< std::complex< double > > > &a, std::vector< int > &ir, int &ifail, std::complex< double > &det, int &jfail)
 
void cfinv (const int n, std::vector< std::vector< std::complex< double > > > &a, std::vector< int > &ir)
 
int cinv (const int n, std::vector< std::vector< std::complex< double > > > &a)
 Replace square matrix A by its inverse.
 

Detailed Description

Linear algebra routines from CERNLIB.

Function Documentation

◆ cfact()

void Garfield::Numerics::CERNLIB::cfact ( const int  n,
std::vector< std::vector< std::complex< double > > > &  a,
std::vector< int > &  ir,
int &  ifail,
std::complex< double > &  det,
int &  jfail 
)

Definition at line 994 of file Numerics.cc.

996 {
997 constexpr double g1 = 1.e-19;
998 constexpr double g2 = 1.e-19;
999
1000 ifail = jfail = 0;
1001
1002 int nxch = 0;
1003 det = std::complex<double>(1., 0.);
1004
1005 for (int j = 1; j <= n; ++j) {
1006 int k = j;
1007 double p = std::max(fabs(real(a[j - 1][j - 1])), fabs(imag(a[j - 1][j - 1])));
1008 if (j == n) {
1009 if (p <= 0.) {
1010 det = std::complex<double>(0., 0.);
1011 ifail = -1;
1012 jfail = 0;
1013 return;
1014 }
1015 det *= a[j - 1][j - 1];
1016 a[j - 1][j - 1] = std::complex<double>(1., 0.) / a[j - 1][j - 1];
1017 const double t = std::max(fabs(real(det)), fabs(imag(det)));
1018 if (t < g1) {
1019 det = std::complex<double>(0., 0.);
1020 if (jfail == 0) jfail = -1;
1021 } else if (t > g2) {
1022 det = std::complex<double>(1., 0.);
1023 if (jfail == 0) jfail = +1;
1024 }
1025 continue;
1026 }
1027 for (int i = j + 1; i <= n; ++i) {
1028 double q = std::max(fabs(real(a[i - 1][j - 1])), fabs(imag(a[i - 1][j - 1])));
1029 if (q <= p) continue;
1030 k = i;
1031 p = q;
1032 }
1033 if (k != j) {
1034 for (int l = 1; l <= n; ++l) {
1035 const auto tf = a[j - 1][l - 1];
1036 a[j - 1][l - 1] = a[k - 1][l - 1];
1037 a[k - 1][l - 1] = tf;
1038 }
1039 ++nxch;
1040 ir[nxch - 1] = j * 4096 + k;
1041 } else if (p <= 0.) {
1042 det = std::complex<double>(0., 0.);
1043 ifail = -1;
1044 jfail = 0;
1045 return;
1046 }
1047 det *= a[j - 1][j - 1];
1048 a[j - 1][j - 1] = 1. / a[j - 1][j - 1];
1049 const double t = std::max(fabs(real(det)), fabs(imag(det)));
1050 if (t < g1) {
1051 det = std::complex<double>(0., 0.);
1052 if (jfail == 0) jfail = -1;
1053 } else if (t > g2) {
1054 det = std::complex<double>(1., 0.);
1055 if (jfail == 0) jfail = +1;
1056 }
1057 for (k = j + 1; k <= n; ++k) {
1058 auto s11 = -a[j - 1][k - 1];
1059 auto s12 = -a[k - 1][j];
1060 if (j == 1) {
1061 a[j - 1][k - 1] = -s11 * a[j - 1][j - 1];
1062 a[k - 1][j] = -(s12 + a[j - 1][j] * a[k - 1][j - 1]);
1063 continue;
1064 }
1065 for (int i = 1; i <= j - 1; ++i) {
1066 s11 += a[i - 1][k - 1] * a[j - 1][i - 1];
1067 s12 += a[i - 1][j] * a[k - 1][i - 1];
1068 }
1069 a[j - 1][k - 1] = -s11 * a[j - 1][j - 1];
1070 a[k - 1][j] = -a[j - 1][j] * a[k - 1][j - 1] - s12;
1071 }
1072 }
1073
1074 if (nxch % 2 != 0) det = -det;
1075 if (jfail != 0) det = std::complex<double>(0., 0.);
1076 ir[n - 1] = nxch;
1077}
Heed::fabs
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:615

Referenced by cinv().

◆ cfinv()

void Garfield::Numerics::CERNLIB::cfinv ( const int  n,
std::vector< std::vector< std::complex< double > > > &  a,
std::vector< int > &  ir 
)

Definition at line 1079 of file Numerics.cc.

1080 {
1081
1082 if (n <= 1) return;
1083 a[1][0] = -a[1][1] * a[0][0] * a[1][0];
1084 a[0][1] = -a[0][1];
1085 if (n > 2) {
1086 for (int i = 3; i <= n; ++i) {
1087 for (int j = 1; j <= i - 2; ++j) {
1088 auto s31 = std::complex<double>(0., 0.);
1089 auto s32 = a[j - 1][i - 1];
1090 for (int k = j; k <= i - 2; ++k) {
1091 s31 += a[k - 1][j - 1] * a[i - 1][k - 1];
1092 s32 += a[j - 1][k] * a[k][i - 1];
1093 }
1094 a[i - 1][j - 1] =
1095 -a[i - 1][i - 1] * (s31 + a[i - 2][j - 1] * a[i - 1][i - 2]);
1096 a[j - 1][i - 1] = -s32;
1097 }
1098 a[i - 1][i - 2] = -a[i - 1][i - 1] * a[i - 2][i - 2] * a[i - 1][i - 2];
1099 a[i - 2][i - 1] = -a[i - 2][i - 1];
1100 }
1101 }
1102
1103 for (int i = 1; i <= n - 1; ++i) {
1104 for (int j = 1; j <= i; ++j) {
1105 auto s33 = a[i - 1][j - 1];
1106 for (int k = 1; k <= n - i; ++k) {
1107 s33 += a[i + k - 1][j - 1] * a[i - 1][i + k - 1];
1108 }
1109 a[i - 1][j - 1] = s33;
1110 }
1111 for (int j = 1; j <= n - i; ++j) {
1112 std::complex<double> s34(0., 0.);
1113 for (int k = j; k <= n - i; ++k) {
1114 s34 += a[i + k - 1][i + j - 1] * a[i - 1][i + k - 1];
1115 }
1116 a[i - 1][i + j - 1] = s34;
1117 }
1118 }
1119
1120 int nxch = ir[n - 1];
1121 if (nxch == 0) return;
1122
1123 for (int m = 1; m <= nxch; ++m) {
1124 int k = nxch - m + 1;
1125 const int ij = ir[k - 1];
1126 const int i = ij / 4096;
1127 const int j = ij % 4096;
1128 for (k = 1; k <= n; ++k) {
1129 std::swap(a[k - 1][i - 1], a[k - 1][j - 1]);
1130 }
1131 }
1132}

Referenced by cinv().

◆ cinv()

int Garfield::Numerics::CERNLIB::cinv ( const int  n,
std::vector< std::vector< std::complex< double > > > &  a 
)

Replace square matrix A by its inverse.

Definition at line 1134 of file Numerics.cc.

1134 {
1135
1136 // Test for parameter errors.
1137 if (n < 1) return 1;
1138
1139 std::complex<double> det(0., 0.);
1140 if (n > 3) {
1141 // n > 3 cases. Factorize matrix and invert.
1142 std::vector<int> ir(n, 0);
1143 int ifail = 0, jfail = 0;
1144 cfact(n, a, ir, ifail, det, jfail);
1145 if (ifail != 0) return ifail;
1146 cfinv(n, a, ir);
1147 } else if (n == 3) {
1148 // n = 3 case. Compute cofactors.
1149 const auto c11 = a[1][1] * a[2][2] - a[1][2] * a[2][1];
1150 const auto c12 = a[1][2] * a[2][0] - a[1][0] * a[2][2];
1151 const auto c13 = a[1][0] * a[2][1] - a[1][1] * a[2][0];
1152 const auto c21 = a[2][1] * a[0][2] - a[2][2] * a[0][1];
1153 const auto c22 = a[2][2] * a[0][0] - a[2][0] * a[0][2];
1154 const auto c23 = a[2][0] * a[0][1] - a[2][1] * a[0][0];
1155 const auto c31 = a[0][1] * a[1][2] - a[0][2] * a[1][1];
1156 const auto c32 = a[0][2] * a[1][0] - a[0][0] * a[1][2];
1157 const auto c33 = a[0][0] * a[1][1] - a[0][1] * a[1][0];
1158 const double t1 = fabs(real(a[0][0])) + fabs(imag(a[0][0]));
1159 const double t2 = fabs(real(a[1][0])) + fabs(imag(a[1][0]));
1160 const double t3 = fabs(real(a[2][0])) + fabs(imag(a[2][0]));
1161
1162 // Set temp = pivot and det = pivot * det.
1163 std::complex<double> temp(0., 0.);
1164 if (t2 < t1 && t3 < t1) {
1165 temp = a[0][0];
1166 det = c22 * c33 - c23 * c32;
1167 } else if (t1 < t2 && t3 < t2) {
1168 temp = a[1][0];
1169 det = c13 * c32 - c12 * c33;
1170 } else {
1171 temp = a[2][0];
1172 det = c23 * c12 - c22 * c13;
1173 }
1174 // Set elements of inverse in A.
1175 if (real(det) == 0. && imag(det) == 0.) return -1;
1176 const auto s = temp / det;
1177 a[0][0] = s * c11;
1178 a[0][1] = s * c21;
1179 a[0][2] = s * c31;
1180 a[1][0] = s * c12;
1181 a[1][1] = s * c22;
1182 a[1][2] = s * c32;
1183 a[2][0] = s * c13;
1184 a[2][1] = s * c23;
1185 a[2][2] = s * c33;
1186 } else if (n == 2) {
1187 // n=2 case by Cramer's rule.
1188 det = a[0][0] * a[1][1] - a[0][1] * a[1][0];
1189 if (real(det) == 0. && imag(det) == 0.) return -1;
1190 const auto s = std::complex<double>(1., 0.) / det;
1191 const auto c11 = s * a[1][1];
1192 a[0][1] = -s * a[0][1];
1193 a[1][0] = -s * a[1][0];
1194 a[1][1] = s * a[0][0];
1195 a[0][0] = c11;
1196 } else {
1197 // n = 1 case.
1198 if (real(a[0][0]) == 0. && imag(a[0][0]) == 0.) return -1;
1199 a[0][0] = std::complex<double>(1., 0.) / a[0][0];
1200 }
1201 return 0;
1202}
Garfield::Numerics::CERNLIB::cfact
void cfact(const int n, std::vector< std::vector< std::complex< double > > > &a, std::vector< int > &ir, int &ifail, std::complex< double > &det, int &jfail)
Definition: Numerics.cc:994
Garfield::Numerics::CERNLIB::cfinv
void cfinv(const int n, std::vector< std::vector< std::complex< double > > > &a, std::vector< int > &ir)
Definition: Numerics.cc:1079

◆ deqinv()

int Garfield::Numerics::CERNLIB::deqinv ( const int  n,
std::vector< std::vector< double > > &  a,
std::vector< double > &  b 
)

Replaces b by the solution x of Ax = b, and replace A by its inverse.

Definition at line 909 of file Numerics.cc.

910 {
911
912 // Test for parameter errors.
913 if (n < 1) return 1;
914
915 double det = 0.;
916 if (n > 3) {
917 // n > 3 cases. Factorize matrix, invert and solve system.
918 std::vector<int> ir(n, 0);
919 int ifail = 0, jfail = 0;
920 dfact(n, a, ir, ifail, det, jfail);
921 if (ifail != 0) return ifail;
922 dfeqn(n, a, ir, b);
923 dfinv(n, a, ir);
924 } else if (n == 3) {
925 // n = 3 case. Compute cofactors.
926 const double c11 = a[1][1] * a[2][2] - a[1][2] * a[2][1];
927 const double c12 = a[1][2] * a[2][0] - a[1][0] * a[2][2];
928 const double c13 = a[1][0] * a[2][1] - a[1][1] * a[2][0];
929 const double c21 = a[2][1] * a[0][2] - a[2][2] * a[0][1];
930 const double c22 = a[2][2] * a[0][0] - a[2][0] * a[0][2];
931 const double c23 = a[2][0] * a[0][1] - a[2][1] * a[0][0];
932 const double c31 = a[0][1] * a[1][2] - a[0][2] * a[1][1];
933 const double c32 = a[0][2] * a[1][0] - a[0][0] * a[1][2];
934 const double c33 = a[0][0] * a[1][1] - a[0][1] * a[1][0];
935 const double t1 = fabs(a[0][0]);
936 const double t2 = fabs(a[1][0]);
937 const double t3 = fabs(a[2][0]);
938
939 // Set temp = pivot and det = pivot * det.
940 double temp = 0.;
941 if (t2 < t1 && t3 < t1) {
942 temp = a[0][0];
943 det = c22 * c33 - c23 * c32;
944 } else if (t1 < t2 && t3 < t2) {
945 temp = a[1][0];
946 det = c13 * c32 - c12 * c33;
947 } else {
948 temp = a[2][0];
949 det = c23 * c12 - c22 * c13;
950 }
951
952 // Set elements of inverse in A.
953 if (det == 0.) return -1;
954 const double s = temp / det;
955 a[0][0] = s * c11;
956 a[0][1] = s * c21;
957 a[0][2] = s * c31;
958 a[1][0] = s * c12;
959 a[1][1] = s * c22;
960 a[1][2] = s * c32;
961 a[2][0] = s * c13;
962 a[2][1] = s * c23;
963 a[2][2] = s * c33;
964
965 // Replace b by Ainv * b.
966 const double b1 = b[0];
967 const double b2 = b[1];
968 b[0] = a[0][0] * b1 + a[0][1] * b2 + a[0][2] * b[2];
969 b[1] = a[1][0] * b1 + a[1][1] * b2 + a[1][2] * b[2];
970 b[2] = a[2][0] * b1 + a[2][1] * b2 + a[2][2] * b[2];
971 } else if (n == 2) {
972 // n = 2 case by Cramer's rule.
973 det = a[0][0] * a[1][1] - a[0][1] * a[1][0];
974 if (det == 0.) return -1;
975 const double s = 1. / det;
976 const double c11 = s * a[1][1];
977 a[0][1] = -s * a[0][1];
978 a[1][0] = -s * a[1][0];
979 a[1][1] = s * a[0][0];
980 a[0][0] = c11;
981
982 const double b1 = b[0];
983 b[0] = c11 * b1 + a[0][1] * b[1];
984 b[1] = a[1][0] * b1 + a[1][1] * b[1];
985 } else {
986 // n = 1 case.
987 if (a[0][0] == 0.) return -1;
988 a[0][0] = 1. / a[0][0];
989 b[0] = a[0][0] * b[0];
990 }
991 return 0;
992}
Garfield::Numerics::CERNLIB::dfinv
void dfinv(const int n, std::vector< std::vector< double > > &a, std::vector< int > &ir)
Definition: Numerics.cc:854
Garfield::Numerics::CERNLIB::dfact
void dfact(const int n, std::vector< std::vector< double > > &a, std::vector< int > &ir, int &ifail, double &det, int &jfail)
Definition: Numerics.cc:668
Garfield::Numerics::CERNLIB::dfeqn
void dfeqn(const int n, std::vector< std::vector< double > > &a, std::vector< int > &ir, std::vector< double > &b)
Definition: Numerics.cc:753

◆ deqn()

int Garfield::Numerics::CERNLIB::deqn ( const int  n,
std::vector< std::vector< double > > &  a,
std::vector< double > &  b 
)

Replaces b by the solution x of Ax = b, after which A is undefined.

Parameters
norder of the square matrix A.
an by n matrix.
bright-hand side vector.
Returns
0: normal exit, -1: singular matrix.

Definition at line 593 of file Numerics.cc.

594 {
595
596 // REPLACES B BY THE SOLUTION X OF A*X=B, AFTER WHICH A IS UNDEFINED.
597
598 if (n < 1) return 1;
599 if (n == 1) {
600 if (a[0][0] == 0.) return -1;
601 const double s = 1. / a[0][0];
602 b[0] *= s;
603 } else if (n == 2) {
604 // Cramer's rule.
605 const double det = a[0][0] * a[1][1] - a[0][1] * a[1][0];
606 if (det == 0.) return -1;
607 const double s = 1. / det;
608 const double b1 = b[0];
609 b[0] = s * ( a[1][1] * b1 - a[0][1] * b[1]);
610 b[1] = s * (-a[1][0] * b1 + a[0][0] * b[1]);
611 } else if (n == 3) {
612 // Factorize matrix A=L*U.
613 // First pivot search.
614 const double t1 = std::abs(a[0][0]);
615 const double t2 = std::abs(a[1][0]);
616 const double t3 = std::abs(a[2][0]);
617 unsigned int m1 = 0, m2 = 0, m3 = 0;
618 if (t1 < t2 && t3 < t2) {
619 // Pivot is A21
620 m1 = 1;
621 m2 = 0;
622 m3 = 2;
623 } else if (t2 < t1 && t3 < t1) {
624 // Pivot is A11
625 m1 = 0;
626 m2 = 1;
627 m3 = 2;
628 } else {
629 // Pivot is A31
630 m1 = 2;
631 m2 = 1;
632 m3 = 0;
633 }
634 double temp = a[m1][0];
635 if (temp == 0.) return deqnGen(n, a, b);
636 const double l11 = 1. / temp;
637 const double u12 = l11 * a[m1][1];
638 const double u13 = l11 * a[m1][2];
639 double l22 = a[m2][1] - a[m2][0] * u12;
640 double l32 = a[m3][1] - a[m3][0] * u12;
641 // Second pivot search.
642 if (std::abs(l22) < std::abs(l32)) {
643 std::swap(m2, m3);
644 std::swap(l22, l32);
645 }
646 double l21 = a[m2][0];
647 double l31 = a[m3][0];
648 if (l22 == 0.) return deqnGen(n, a, b);
649 l22 = 1. / l22;
650 const double u23 = l22 * (a[m2][2] - l21 * u13);
651 temp = a[m3][2] - l31 * u13 - l32 * u23;
652 if (temp == 0.) return deqnGen(n, a, b);
653 const double l33 = 1. / temp;
654
655 // Solve L*Y=B and U*X=Y.
656 const double y1 = l11 * b[m1];
657 const double y2 = l22 * (b[m2] - l21 * y1);
658 b[2] = l33 * (b[m3] - l31 * y1 - l32 * y2);
659 b[1] = y2 - u23 * b[2];
660 b[0] = y1 - u12 * b[1] - u13 * b[2];
661 } else {
662 // N > 3 cases. Factorize matrix and solve system.
663 return deqnGen(n, a, b);
664 }
665 return 0;
666}

◆ dfact()

void Garfield::Numerics::CERNLIB::dfact ( const int  n,
std::vector< std::vector< double > > &  a,
std::vector< int > &  ir,
int &  ifail,
double &  det,
int &  jfail 
)

Definition at line 668 of file Numerics.cc.

669 {
670 constexpr double g1 = 1.e-19;
671 constexpr double g2 = 1.e-19;
672
673 double t;
674 int k;
675
676 ifail = jfail = 0;
677
678 int nxch = 0;
679 det = 1.;
680
681 for (int j = 1; j <= n; ++j) {
682 k = j;
683 double p = fabs(a[j - 1][j - 1]);
684 if (j == n) {
685 if (p <= 0.) {
686 det = 0.;
687 ifail = -1;
688 jfail = 0;
689 return;
690 }
691 det *= a[j - 1][j - 1];
692 a[j - 1][j - 1] = 1. / a[j - 1][j - 1];
693 t = fabs(det);
694 if (t < g1) {
695 det = 0.;
696 if (jfail == 0) jfail = -1;
697 } else if (t > g2) {
698 det = 1.;
699 if (jfail == 0) jfail = +1;
700 }
701 continue;
702 }
703 for (int i = j + 1; i <= n; ++i) {
704 double q = std::abs(a[i - 1][j - 1]);
705 if (q <= p) continue;
706 k = i;
707 p = q;
708 }
709 if (k != j) {
710 for (int l = 1; l <= n; ++l) {
711 std::swap(a[j - 1][l - 1], a[k - 1][l - 1]);
712 }
713 ++nxch;
714 ir[nxch - 1] = j * 4096 + k;
715 } else if (p <= 0.) {
716 det = 0.;
717 ifail = -1;
718 jfail = 0;
719 return;
720 }
721 det *= a[j - 1][j - 1];
722 a[j - 1][j - 1] = 1. / a[j - 1][j - 1];
723 t = fabs(det);
724 if (t < g1) {
725 det = 0.;
726 if (jfail == 0) jfail = -1;
727 } else if (t > g2) {
728 det = 1.;
729 if (jfail == 0) jfail = +1;
730 }
731 for (k = j + 1; k <= n; ++k) {
732 double s11 = -a[j - 1][k - 1];
733 double s12 = -a[k - 1][j];
734 if (j == 1) {
735 a[j - 1][k - 1] = -s11 * a[j - 1][j - 1];
736 a[k - 1][j] = -(s12 + a[j - 1][j] * a[k - 1][j - 1]);
737 continue;
738 }
739 for (int i = 1; i <= j - 1; ++i) {
740 s11 += a[i - 1][k - 1] * a[j - 1][i - 1];
741 s12 += a[i - 1][j] * a[k - 1][i - 1];
742 }
743 a[j - 1][k - 1] = -s11 * a[j - 1][j - 1];
744 a[k - 1][j] = -a[j - 1][j] * a[k - 1][j - 1] - s12;
745 }
746 }
747
748 if (nxch % 2 != 0) det = -det;
749 if (jfail != 0) det = 0.;
750 ir[n - 1] = nxch;
751}

Referenced by deqinv(), and dinv().

◆ dfeqn()

void Garfield::Numerics::CERNLIB::dfeqn ( const int  n,
std::vector< std::vector< double > > &  a,
std::vector< int > &  ir,
std::vector< double > &  b 
)

Definition at line 753 of file Numerics.cc.

754 {
755 if (n <= 0) return;
756
757 int nxch = ir[n - 1];
758 if (nxch != 0) {
759 for (int m = 1; m <= nxch; ++m) {
760 const int ij = ir[m - 1];
761 const int i = ij / 4096;
762 const int j = ij % 4096;
763 std::swap(b[i - 1], b[j - 1]);
764 }
765 }
766
767 b[0] *= a[0][0];
768 if (n == 1) return;
769
770 for (int i = 2; i <= n; ++i) {
771 double s21 = -b[i - 1];
772 for (int j = 1; j <= i - 1; ++j) {
773 s21 += a[i - 1][j - 1] * b[j - 1];
774 }
775 b[i - 1] = -a[i - 1][i - 1] * s21;
776 }
777
778 for (int i = 1; i <= n - 1; ++i) {
779 double s22 = -b[n - i - 1];
780 for (int j = 1; j <= i; ++j) {
781 s22 += a[n - i - 1][n - j] * b[n - j];
782 }
783 b[n - i - 1] = -s22;
784 }
785}

Referenced by deqinv().

◆ dfinv()

void Garfield::Numerics::CERNLIB::dfinv ( const int  n,
std::vector< std::vector< double > > &  a,
std::vector< int > &  ir 
)

Definition at line 854 of file Numerics.cc.

855 {
856
857 if (n <= 1) return;
858 a[1][0] = -a[1][1] * a[0][0] * a[1][0];
859 a[0][1] = -a[0][1];
860 if (n > 2) {
861 for (int i = 3; i <= n; ++i) {
862 for (int j = 1; j <= i - 2; ++j) {
863 double s31 = 0.;
864 double s32 = a[j - 1][i - 1];
865 for (int k = j; k <= i - 2; ++k) {
866 s31 += a[k - 1][j - 1] * a[i - 1][k - 1];
867 s32 += a[j - 1][k] * a[k][i - 1];
868 }
869 a[i - 1][j - 1] =
870 -a[i - 1][i - 1] * (s31 + a[i - 2][j - 1] * a[i - 1][i - 2]);
871 a[j - 1][i - 1] = -s32;
872 }
873 a[i - 1][i - 2] = -a[i - 1][i - 1] * a[i - 2][i - 2] * a[i - 1][i - 2];
874 a[i - 2][i - 1] = -a[i - 2][i - 1];
875 }
876 }
877
878 for (int i = 1; i <= n - 1; ++i) {
879 for (int j = 1; j <= i; ++j) {
880 double s33 = a[i - 1][j - 1];
881 for (int k = 1; k <= n - i; ++k) {
882 s33 += a[i + k - 1][j - 1] * a[i - 1][i + k - 1];
883 }
884 a[i - 1][j - 1] = s33;
885 }
886 for (int j = 1; j <= n - i; ++j) {
887 double s34 = 0.;
888 for (int k = j; k <= n - i; ++k) {
889 s34 += a[i + k - 1][i + j - 1] * a[i - 1][i + k - 1];
890 }
891 a[i - 1][i + j - 1] = s34;
892 }
893 }
894
895 int nxch = ir[n - 1];
896 if (nxch == 0) return;
897
898 for (int m = 1; m <= nxch; ++m) {
899 int k = nxch - m + 1;
900 int ij = ir[k - 1];
901 int i = ij / 4096;
902 int j = ij % 4096;
903 for (k = 1; k <= n; ++k) {
904 std::swap(a[k - 1][i - 1], a[k - 1][j - 1]);
905 }
906 }
907}

Referenced by deqinv(), and dinv().

◆ dinv()

int Garfield::Numerics::CERNLIB::dinv ( const int  n,
std::vector< std::vector< double > > &  a 
)

Replace square matrix A by its inverse.

Definition at line 787 of file Numerics.cc.

787 {
788
789 if (n < 1) return 1;
790 if (n > 3) {
791 // Factorize matrix and invert.
792 double det = 0.;
793 int ifail = 0;
794 int jfail = 0;
795 std::vector<int> ir(n, 0);
796 dfact(n, a, ir, ifail, det, jfail);
797 if (ifail != 0) return ifail;
798 dfinv(n, a, ir);
799 } else if (n == 3) {
800 // Compute cofactors.
801 const double c11 = a[1][1] * a[2][2] - a[1][2] * a[2][1];
802 const double c12 = a[1][2] * a[2][0] - a[1][0] * a[2][2];
803 const double c13 = a[1][0] * a[2][1] - a[1][1] * a[2][0];
804 const double c21 = a[2][1] * a[0][2] - a[2][2] * a[0][1];
805 const double c22 = a[2][2] * a[0][0] - a[2][0] * a[0][2];
806 const double c23 = a[2][0] * a[0][1] - a[2][1] * a[0][0];
807 const double c31 = a[0][1] * a[1][2] - a[0][2] * a[1][1];
808 const double c32 = a[0][2] * a[1][0] - a[0][0] * a[1][2];
809 const double c33 = a[0][0] * a[1][1] - a[0][1] * a[1][0];
810 const double t1 = fabs(a[0][0]);
811 const double t2 = fabs(a[1][0]);
812 const double t3 = fabs(a[2][0]);
813 double det = 0.;
814 double pivot = 0.;
815 if (t2 < t1 && t3 < t1) {
816 pivot = a[0][0];
817 det = c22 * c33 - c23 * c32;
818 } else if (t1 < t2 && t3 < t2) {
819 pivot = a[1][0];
820 det = c13 * c32 - c12 * c33;
821 } else {
822 pivot = a[2][0];
823 det = c23 * c12 - c22 * c13;
824 }
825 // Set elements of inverse in A.
826 if (det == 0.) return -1;
827 double s = pivot / det;
828 a[0][0] = s * c11;
829 a[0][1] = s * c21;
830 a[0][2] = s * c31;
831 a[1][0] = s * c12;
832 a[1][1] = s * c22;
833 a[1][2] = s * c32;
834 a[2][0] = s * c13;
835 a[2][1] = s * c23;
836 a[2][2] = s * c33;
837 } else if (n == 2) {
838 // Cramer's rule.
839 const double det = a[0][0] * a[1][1] - a[0][1] * a[1][0];
840 if (det == 0.) return -1;
841 const double s = 1. / det;
842 const double c11 = s * a[1][1];
843 a[0][1] = -s * a[0][1];
844 a[1][0] = -s * a[1][0];
845 a[1][1] = s * a[0][0];
846 a[0][0] = c11;
847 } else if (n == 1) {
848 if (a[0][0] == 0.) return -1;
849 a[0][0] = 1. / a[0][0];
850 }
851 return 0;
852}

Referenced by Garfield::Numerics::LeastSquaresFit().