whitespace
This commit is contained in:
@ -1247,7 +1247,7 @@ void FitPOD::least_squares_fit(datastruct data)
|
|||||||
if (comm->me == 0) {
|
if (comm->me == 0) {
|
||||||
|
|
||||||
// symmetrize A
|
// symmetrize A
|
||||||
|
|
||||||
for (int i = 0; i<nd; i++)
|
for (int i = 0; i<nd; i++)
|
||||||
for (int j = i; j<nd; j++) {
|
for (int j = i; j<nd; j++) {
|
||||||
double a1 = desc.A[i + nd*j];
|
double a1 = desc.A[i + nd*j];
|
||||||
@ -1255,20 +1255,20 @@ void FitPOD::least_squares_fit(datastruct data)
|
|||||||
desc.A[i + nd*j] = 0.5*(a1+a2);
|
desc.A[i + nd*j] = 0.5*(a1+a2);
|
||||||
desc.A[j + nd*i] = 0.5*(a1+a2);
|
desc.A[j + nd*i] = 0.5*(a1+a2);
|
||||||
}
|
}
|
||||||
|
|
||||||
// scale A and b
|
// scale A and b
|
||||||
|
|
||||||
double maxb = 0.0;
|
double maxb = 0.0;
|
||||||
for (int i = 0; i<nd; i++)
|
for (int i = 0; i<nd; i++)
|
||||||
maxb = (maxb > fabs(desc.b[i])) ? maxb : fabs(desc.b[i]);
|
maxb = (maxb > fabs(desc.b[i])) ? maxb : fabs(desc.b[i]);
|
||||||
|
|
||||||
maxb = 1.0/maxb;
|
maxb = 1.0/maxb;
|
||||||
for (int i = 0; i<nd; i++)
|
for (int i = 0; i<nd; i++)
|
||||||
desc.b[i] = desc.b[i]*maxb;
|
desc.b[i] = desc.b[i]*maxb;
|
||||||
|
|
||||||
for (int i = 0; i<nd*nd; i++)
|
for (int i = 0; i<nd*nd; i++)
|
||||||
desc.A[i] = desc.A[i]*maxb;
|
desc.A[i] = desc.A[i]*maxb;
|
||||||
|
|
||||||
for (int i = 0; i<nd; i++) {
|
for (int i = 0; i<nd; i++) {
|
||||||
desc.c[i] = desc.b[i];
|
desc.c[i] = desc.b[i];
|
||||||
desc.A[i + nd*i] = desc.A[i + nd*i]*(1.0 + SMALL);
|
desc.A[i + nd*i] = desc.A[i + nd*i]*(1.0 + SMALL);
|
||||||
@ -1276,13 +1276,13 @@ void FitPOD::least_squares_fit(datastruct data)
|
|||||||
}
|
}
|
||||||
|
|
||||||
// solving the linear system A * c = b
|
// solving the linear system A * c = b
|
||||||
|
|
||||||
int nrhs=1, info;
|
int nrhs=1, info;
|
||||||
char chu = 'U';
|
char chu = 'U';
|
||||||
DPOSV(&chu, &nd, &nrhs, desc.A, &nd, desc.c, &nd, &info);
|
DPOSV(&chu, &nd, &nrhs, desc.A, &nd, desc.c, &nd, &info);
|
||||||
|
|
||||||
// compute the inverse of A
|
// compute the inverse of A
|
||||||
|
|
||||||
// int lwork = 10*nd;
|
// int lwork = 10*nd;
|
||||||
// int info;
|
// int info;
|
||||||
// int *ipiv;
|
// int *ipiv;
|
||||||
@ -1293,14 +1293,14 @@ void FitPOD::least_squares_fit(datastruct data)
|
|||||||
// DGETRI(&nd, desc.A, &nd, ipiv, work, &lwork, &info);
|
// DGETRI(&nd, desc.A, &nd, ipiv, work, &lwork, &info);
|
||||||
// memory->destroy(ipiv);
|
// memory->destroy(ipiv);
|
||||||
// memory->destroy(work);
|
// memory->destroy(work);
|
||||||
|
|
||||||
// compute c = inverse(A) * b
|
// compute c = inverse(A) * b
|
||||||
|
|
||||||
// for (int i = 0; i<nd; i++) {
|
// for (int i = 0; i<nd; i++) {
|
||||||
// desc.c[i] = 0.0;
|
// desc.c[i] = 0.0;
|
||||||
// for (int j = 0; j<nd; j++)
|
// for (int j = 0; j<nd; j++)
|
||||||
// desc.c[i] += desc.A[i + nd*j]*desc.b[j];
|
// desc.c[i] += desc.A[i + nd*j]*desc.b[j];
|
||||||
// }
|
// }
|
||||||
}
|
}
|
||||||
|
|
||||||
MPI_Bcast(desc.c, nd, MPI_DOUBLE, 0, world);
|
MPI_Bcast(desc.c, nd, MPI_DOUBLE, 0, world);
|
||||||
@ -1311,13 +1311,13 @@ void FitPOD::least_squares_fit(datastruct data)
|
|||||||
std::string filename = podptr->pod.filenametag + "_coefficients" + ".pod";
|
std::string filename = podptr->pod.filenametag + "_coefficients" + ".pod";
|
||||||
FILE *fp = fopen(filename.c_str(), "w");
|
FILE *fp = fopen(filename.c_str(), "w");
|
||||||
|
|
||||||
//int prec = utils::inumeric(FLERR,podptr->pod.precision,false,lmp);
|
//int prec = utils::inumeric(FLERR,podptr->pod.precision,false,lmp);
|
||||||
std::string prec = "{:<10." + podptr->pod.precision + "f}\n";
|
std::string prec = "{:<10." + podptr->pod.precision + "f}\n";
|
||||||
|
|
||||||
fmt::print(fp, "POD_coefficients: {}\n", nd);
|
fmt::print(fp, "POD_coefficients: {}\n", nd);
|
||||||
for (int count = 0; count < nd; count++) {
|
for (int count = 0; count < nd; count++) {
|
||||||
fmt::print(fp, prec.c_str(), desc.c[count]);
|
fmt::print(fp, prec.c_str(), desc.c[count]);
|
||||||
}
|
}
|
||||||
fclose(fp);
|
fclose(fp);
|
||||||
utils::logmesg(lmp, "**************** End of Least-Squares Fitting ****************\n");
|
utils::logmesg(lmp, "**************** End of Least-Squares Fitting ****************\n");
|
||||||
}
|
}
|
||||||
|
|||||||
Reference in New Issue
Block a user