Garfield++ 5.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.
 
void cfft (std::vector< std::complex< double > > &a, const int msign)
 

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 995 of file Numerics.cc.

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

Referenced by cinv().

◆ cfft()

void Garfield::Numerics::CERNLIB::cfft ( std::vector< std::complex< double > > & a,
const int msign )

Definition at line 1205 of file Numerics.cc.

1205 {
1206
1207 if (msign == 0) return;
1208 const int m = std::abs(msign);
1209 const int n = pow(2, m);
1210 // Bit reversal.
1211 int j = 1;
1212 for (int i = 1; i <= n - 1; ++i) {
1213 if (i < j) std::swap(a[j - 1], a[i - 1]);
1214 int k = n / 2;
1215 while (k < j) {
1216 j -= k;
1217 k /= 2;
1218 }
1219 j += k;
1220 }
1221 for (int i = 1; i <= n; i += 2) {
1222 const auto t = a[i];
1223 a[i] = a[i - 1] - t;
1224 a[i - 1] += t;
1225 }
1226 if (m == 1) return;
1227 double c = 0.;
1228 double s = msign >= 0 ? 1. : -1.;
1229 int le = 2;
1230 for (int l = 2; l <= m; ++l) {
1231 std::complex<double> w(c, s);
1232 std::complex<double> u = w;
1233 c = sqrt(0.5 * c + 0.5);
1234 s = imag(w) / (c + c);
1235 int le1 = le;
1236 le = le1 + le1;
1237 for (int i = 1; i <= n; i += le) {
1238 const auto t = a[i + le1 - 1];
1239 a[i + le1 - 1] = a[i - 1] - t;
1240 a[i - 1] += t;
1241 }
1242 for (j = 2; j <= le1; ++j) {
1243 for (int i = j; i <= n; i += le) {
1244 const auto t = a[i + le1 - 1] * u;
1245 a[i + le1 - 1] = a[i - 1] - t;
1246 a[i - 1] += t;
1247 }
1248 u *= w;
1249 }
1250 }
1251}

◆ cfinv()

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

Definition at line 1080 of file Numerics.cc.

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

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 1135 of file Numerics.cc.

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

◆ 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 910 of file Numerics.cc.

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

◆ 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 594 of file Numerics.cc.

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

◆ 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 669 of file Numerics.cc.

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

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 754 of file Numerics.cc.

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

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 855 of file Numerics.cc.

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

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 788 of file Numerics.cc.

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

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