982 static const int max_array = 25;
984 static std::vector<double> xvec (max_array);
985 static std::vector<int> pivv (max_array);
986 typedef std::vector<int>::iterator pivIter;
988 static std::vector<double,Alloc<double,25> > xvec (max_array);
989 static std::vector<int, Alloc<int, 25> > pivv (max_array);
990 typedef std::vector<int,Alloc<int,25> >::iterator pivIter;
992 if (xvec.size() <
static_cast<unsigned int>(nrow)) xvec.resize(nrow);
993 if (pivv.size() <
static_cast<unsigned int>(nrow)) pivv.resize(nrow);
998 mIter x = xvec.begin();
1000 pivIter piv = pivv.begin();
1003 double temp1, temp2;
1005 double lambda, sigma;
1006 const double alpha = .6404;
1007 const double epsilon = 32*DBL_EPSILON;
1013 for (i = 0; i < nrow; ++i) piv[i] = i+1;
1021 for (j=1; j < nrow; j+=is)
1023 mjj = m.begin() + j*(j-1)/2 + j-1;
1027 for (i=j+1; i <= nrow ; ++i) {
1029 ip = m.begin() + (i-1)*i/2 + j-1;
1030 if (fabs(*ip) > lambda) {
1043 if (fabs(*mjj) >= lambda*alpha) {
1048 ip = m.begin() + pivrow*(pivrow-1)/2+j-1;
1049 for (k=j; k < pivrow; k++) {
1050 if (fabs(*ip) > sigma) sigma = fabs(*ip);
1053 if (sigma * fabs(*mjj) >= alpha * lambda * lambda) {
1056 }
else if (fabs(*(m.begin()+pivrow*(pivrow-1)/2+pivrow-1))
1069 temp2 = *mjj = 1./ *mjj;
1072 for (i=j+1; i <= nrow; i++) {
1073 temp1 = *(m.begin() + i*(i-1)/2 + j-1) * temp2;
1074 ip = m.begin()+i*(i-1)/2+j;
1075 for (k=j+1; k<=i; k++) {
1076 *ip -= temp1 * *(m.begin() + k*(k-1)/2 + j-1);
1077 if (fabs(*ip) <= epsilon)
1084 for (i=j+1; i <= nrow; ++i) {
1086 ip = m.begin() + (i-1)*i/2 + j-1;
1094 ip = m.begin() + pivrow*(pivrow-1)/2 + j;
1095 for (i=j+1; i < pivrow; i++, ip++) {
1096 temp1 = *(m.begin() + i*(i-1)/2 + j-1);
1097 *(m.begin() + i*(i-1)/2 + j-1)= *ip;
1101 *mjj = *(m.begin()+pivrow*(pivrow-1)/2+pivrow-1);
1102 *(m.begin()+pivrow*(pivrow-1)/2+pivrow-1) = temp1;
1103 ip = m.begin() + (pivrow+1)*pivrow/2 + j-1;
1105 for (i = pivrow+1; i <= nrow; ip += i, iq += i++) {
1115 temp2 = *mjj = 1./ *mjj;
1118 for (i = j+1; i <= nrow; i++) {
1119 temp1 = *(m.begin() + i*(i-1)/2 + j-1) * temp2;
1120 ip = m.begin()+i*(i-1)/2+j;
1121 for (k=j+1; k<=i; k++) {
1122 *ip -= temp1 * *(m.begin() + k*(k-1)/2 + j-1);
1123 if (fabs(*ip) <= epsilon)
1130 for (i=j+1; i <= nrow; ++i) {
1132 ip = m.begin() + (i-1)*i/2 + j-1;
1139 if (j+1 != pivrow) {
1142 ip = m.begin() + pivrow*(pivrow-1)/2 + j+1;
1143 for (i=j+2; i < pivrow; i++, ip++) {
1144 temp1 = *(m.begin() + i*(i-1)/2 + j);
1145 *(m.begin() + i*(i-1)/2 + j) = *ip;
1148 temp1 = *(mjj + j + 1);
1150 *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1);
1151 *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1) = temp1;
1153 *(mjj + j) = *(m.begin() + pivrow*(pivrow-1)/2 + j-1);
1154 *(m.begin() + pivrow*(pivrow-1)/2 + j-1) = temp1;
1155 ip = m.begin() + (pivrow+1)*pivrow/2 + j;
1156 iq = ip + pivrow-(j+1);
1157 for (i = pivrow+1; i <= nrow; ip += i, iq += i++) {
1164 temp2 = *mjj * *(mjj + j + 1) - *(mjj + j) * *(mjj + j);
1167 <<
"SymMatrix::bunch_invert: error in pivot choice"
1174 *mjj = *(mjj + j + 1) * temp2;
1175 *(mjj + j + 1) = temp1 * temp2;
1176 *(mjj + j) = - *(mjj + j) * temp2;
1180 for (i=j+2; i <= nrow ; i++) {
1181 ip = m.begin() + i*(i-1)/2 + j-1;
1182 temp1 = *ip * *mjj + *(ip + 1) * *(mjj + j);
1183 if (fabs(temp1 ) <= epsilon) temp1 = 0;
1184 temp2 = *ip * *(mjj + j) + *(ip + 1) * *(mjj + j + 1);
1185 if (fabs(temp2 ) <= epsilon) temp2 = 0;
1186 for (k = j+2; k <= i ; k++) {
1187 ip = m.begin() + i*(i-1)/2 + k-1;
1188 iq = m.begin() + k*(k-1)/2 + j-1;
1189 *ip -= temp1 * *iq + temp2 * *(iq+1);
1190 if (fabs(*ip) <= epsilon)
1195 for (i=j+2; i <= nrow ; i++) {
1196 ip = m.begin() + i*(i-1)/2 + j-1;
1197 temp1 = *ip * *mjj + *(ip+1) * *(mjj + j);
1198 if (fabs(temp1) <= epsilon) temp1 = 0;
1199 *(ip+1) = *ip * *(mjj + j)
1200 + *(ip+1) * *(mjj + j + 1);
1201 if (fabs(*(ip+1)) <= epsilon) *(ip+1) = 0;
1210 mjj = m.begin() + j*(j-1)/2 + j-1;
1214 }
else { *mjj = 1. / *mjj; }
1219 for (j = nrow ; j >= 1 ; j -= is)
1221 mjj = m.begin() + j*(j-1)/2 + j-1;
1227 ip = m.begin() + (j+1)*j/2 - 1;
1228 for (i=0; i < nrow-j; ++i) {
1232 for (i=j+1; i<=nrow ; i++) {
1234 ip = m.begin() + i*(i-1)/2 + j;
1235 for (k=0; k <= i-j-1; k++) temp2 += *ip++ * x[k];
1239 for ( ; k < nrow-j; ++k) {
1241 temp2 += *ip * x[k];
1243 *(m.begin()+ i*(i-1)/2 + j-1) = -temp2;
1249 ip = m.begin() + (j+1)*j/2 - 1;
1250 for (k=0; k < nrow-j; ++k) {
1252 temp2 += x[k] * *ip;
1258 std::cerr <<
"error in piv" << piv[j-1] << std::endl;
1263 ip = m.begin() + (j+1)*j/2 - 1;
1264 for (i=0; i < nrow-j; ++i) {
1268 for (i=j+1; i<=nrow ; i++) {
1270 ip = m.begin() + i*(i-1)/2 + j;
1271 for (k=0; k <= i-j-1; k++)
1272 temp2 += *ip++ * x[k];
1273 for (ip += i-1; k < nrow-j; ip += 1+j+k++)
1274 temp2 += *ip * x[k];
1275 *(m.begin()+ i*(i-1)/2 + j-1) = -temp2;
1280 ip = m.begin() + (j+1)*j/2 - 1;
1281 for (k=0; k < nrow-j; ++k) {
1283 temp2 += x[k] * *ip;
1289 ip = m.begin() + (j+1)*j/2 - 2;
1290 for (i=j+1; i <= nrow; ++i) {
1292 temp2 += *ip * *(ip+1);
1297 ip = m.begin() + (j+1)*j/2 - 2;
1298 for (i=0; i < nrow-j; ++i) {
1302 for (i=j+1; i <= nrow ; i++) {
1304 ip = m.begin() + i*(i-1)/2 + j;
1305 for (k=0; k <= i-j-1; k++)
1306 temp2 += *ip++ * x[k];
1307 for (ip += i-1; k < nrow-j; ip += 1+j+k++)
1308 temp2 += *ip * x[k];
1309 *(m.begin()+ i*(i-1)/2 + j-2)= -temp2;
1315 ip = m.begin() + (j+1)*j/2 - 2;
1316 for (k=0; k < nrow-j; ++k) {
1318 temp2 += x[k] * *ip;
1327 pivrow = (piv[j-1]==0)? -piv[j-2] : piv[j-1];
1328 ip = m.begin() + pivrow*(pivrow-1)/2 + j;
1329 for (i=j+1;i < pivrow; i++, ip++) {
1330 temp1 = *(m.begin() + i*(i-1)/2 + j-1);
1331 *(m.begin() + i*(i-1)/2 + j-1) = *ip;
1335 *mjj = *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1);
1336 *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1) = temp1;
1339 *(mjj-1) = *( m.begin() + pivrow*(pivrow-1)/2 + j-2);
1340 *( m.begin() + pivrow*(pivrow-1)/2 + j-2) = temp1;
1344 if( pivrow < nrow ) {
1345 ip = m.begin() + (pivrow+1)*pivrow/2 + j-1;
1347 iq = ip + (pivrow-j);
1348 for (i = pivrow+1; i <= nrow; i++) {