From d33f13003742f45688a9a481cd7b4814b5f24ac9 Mon Sep 17 00:00:00 2001 From: "joaquin.f.fernandez" Date: Wed, 26 Apr 2023 10:25:02 -0300 Subject: [PATCH] [iss-221] Squashed commit of the following: commit 114b09a3749dae03b87933781bea1463975323e3 Author: joaquin.f.fernandez Date: Tue Apr 25 14:12:11 2023 -0300 Move QSS_TEST method to the end of the Solver enum. commit 1b50762fd09a4a40c1224e8bf1b9cc9fd57eba98 Author: joaquin.f.fernandez Date: Thu Apr 20 10:46:13 2023 -0300 Add qss test method to quantizer. commit 3323d731d9d352b9fb453a882635704a535167b6 Author: joaquin.f.fernandez Date: Thu Apr 20 10:45:57 2023 -0300 Add qss test method. commit 34d36d6d4a090a6a9688395b85ef8f8c498cbe98 Author: joaquin.f.fernandez Date: Thu Apr 20 10:45:32 2023 -0300 Parse qss test methods in settings. commit 5395fcfdf1f631940b7ea0c5eca6efa8ebb18f9d Author: joaquin.f.fernandez Date: Thu Apr 20 10:45:15 2023 -0300 Add qss test methods to data enum. --- src/engine/common/data.h | 3 +- src/engine/common/settings.c | 2 + src/engine/qss/methods/Makefile.include | 3 +- src/engine/qss/methods/qss_test.c | 487 ++++++++++++++++++ src/engine/qss/methods/qss_test.h | 76 +++ .../qss/parallel/methods/Makefile.include | 3 +- src/engine/qss/qss_quantizer.c | 18 + src/engine/qss/qss_quantizer.h | 192 ++----- 8 files changed, 639 insertions(+), 145 deletions(-) create mode 100644 src/engine/qss/methods/qss_test.c create mode 100644 src/engine/qss/methods/qss_test.h diff --git a/src/engine/common/data.h b/src/engine/common/data.h index d4403e8d..c3f244e6 100644 --- a/src/engine/common/data.h +++ b/src/engine/common/data.h @@ -61,7 +61,8 @@ typedef enum { SD_DOPRI, SD_CVODE_BDF, SD_CVODE_AM, - SD_IDA + SD_IDA, + SD_QSS_TEST } SD_Solver; typedef enum { diff --git a/src/engine/common/settings.c b/src/engine/common/settings.c index 102eec34..d6a06257 100644 --- a/src/engine/common/settings.c +++ b/src/engine/common/settings.c @@ -46,6 +46,8 @@ SD_Solver _getSolver(const char *sol) return SD_LIQSS3; } else if (!strcmp(sol, "QSS4")) { return SD_QSS4; + } else if (!strcmp(sol, "QSS_TEST")) { + return SD_QSS_TEST; } else if (!strcmp(sol, "DASSL")) { return SD_DASSL; } else if (!strcmp(sol, "DOPRI")) { diff --git a/src/engine/qss/methods/Makefile.include b/src/engine/qss/methods/Makefile.include index 52094daf..5fba8edf 100644 --- a/src/engine/qss/methods/Makefile.include +++ b/src/engine/qss/methods/Makefile.include @@ -10,7 +10,8 @@ SEQ_METHODS_SRC = $(QSS_METHODS_DIR)/qss.c \ $(QSS_METHODS_DIR)/liqss2.c \ $(QSS_METHODS_DIR)/liqss_bdf.c \ $(QSS_METHODS_DIR)/liqss3.c \ - $(QSS_METHODS_DIR)/cqss.c + $(QSS_METHODS_DIR)/cqss.c \ + $(QSS_METHODS_DIR)/qss_test.c # Objects diff --git a/src/engine/qss/methods/qss_test.c b/src/engine/qss/methods/qss_test.c new file mode 100644 index 00000000..6ed6728b --- /dev/null +++ b/src/engine/qss/methods/qss_test.c @@ -0,0 +1,487 @@ +/***************************************************************************** + + This file is part of QSS Solver. + + QSS Solver is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + QSS Solver is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with QSS Solver. If not, see . + +******************************************************************************/ + +#include +#include +#include + +#include +#include + +#define TOL 2 + +#ifdef QSS_PARALLEL +void QSS_TEST_PAR_init(QA_quantizer quantizer, QSS_data simData, QSS_time simTime) +#else +void QSS_TEST_init(QA_quantizer quantizer, QSS_data simData, QSS_time simTime) +#endif +{ + int states = simData->states; + int i, j; + quantizer->state->oldDx = (double *)malloc(states * sizeof(double)); + quantizer->state->qAux = (double *)malloc(states * sizeof(double)); + quantizer->state->ltq = (double *)malloc(states * sizeof(double)); + quantizer->state->A = (double **)malloc(states * sizeof(double *)); + quantizer->state->U0 = (double **)malloc(states * sizeof(double *)); + quantizer->state->U1 = (double **)malloc(states * sizeof(double *)); + for (i = 0; i < states; i++) { + quantizer->state->A[i] = (double *)malloc(states * sizeof(double)); + quantizer->state->U0[i] = (double *)malloc(states * sizeof(double)); + quantizer->state->U1[i] = (double *)malloc(states * sizeof(double)); + } + quantizer->state->change = (bool *)malloc(states * sizeof(bool)); + quantizer->state->next = (double *)malloc(3 * states * sizeof(double)); + quantizer->state->tx = (double *)malloc(states * sizeof(double)); + + quantizer->state->cont = (double *)malloc(10 * sizeof(double)); + for (i = 0; i < 10; i++) quantizer->state->cont[i] = 0; + + quantizer->state->finTime = simData->ft; + for (i = 0; i < states; i++) { + int cf0 = i * 3; + simData->x[cf0 + 2] = 0; + simData->q[cf0] = simData->x[cf0]; + simData->q[cf0 + 1] = 0; + quantizer->state->qAux[i] = simData->x[cf0]; + quantizer->state->oldDx[i] = 0; + quantizer->state->ltq[i] = simData->it; + quantizer->state->change[i] = FALSE; + quantizer->state->tx[i] = 0; + quantizer->state->next[cf0] = 0; + quantizer->state->next[cf0 + 1] = 0; + quantizer->state->next[cf0 + 2] = 0; + for (j = 0; j < states; j++) { + quantizer->state->A[i][j] = 0; + quantizer->state->U0[i][j] = 0; + quantizer->state->U1[i][j] = 0; + } + } +#ifdef QSS_PARALLEL + quantizer->state->qMap = simData->lp->qMap; +#endif + quantizer->state->minStep = simData->params->minStep; + quantizer->state->lSimTime = simTime; + quantizer->state->nSD = simData->nSD; + quantizer->state->SD = simData->SD; + quantizer->ops->recomputeNextTimes = QSS_TEST_recomputeNextTimes; + quantizer->ops->recomputeNextTime = QSS_TEST_recomputeNextTime; + quantizer->ops->nextTime = QSS_TEST_nextTime; + quantizer->ops->updateQuantizedState = QSS_TEST_updateQuantizedState; +} + +#ifdef QSS_PARALLEL +void QSS_TEST_PAR_recomputeNextTime(QA_quantizer quantizer, int i, double t, double *nTime, double *x, double *lqu, double *q) +#else +void QSS_TEST_recomputeNextTime(QA_quantizer quantizer, int i, double t, double *nTime, double *x, double *lqu, double *q) +#endif +{ + int i0 = i * 3, i1 = i0 + 1, i2 = i1 + 1; + double coeff[3], diffQ, timeaux; + double **A = quantizer->state->A; + double **U0 = quantizer->state->U0; + double **U1 = quantizer->state->U1; + int stind = quantizer->state->lSimTime->minIndex; + bool *change = quantizer->state->change; + bool self = quantizer->state->lSimTime->minIndex == i && quantizer->state->lSimTime->type == ST_State; + double *tx = quantizer->state->tx; + // double *c = quantizer->state->cont; + + if (t > 0) { + diffQ = q[3 * stind] - quantizer->state->qAux[stind]; + if (diffQ) + A[i][stind] = (x[i1] - quantizer->state->oldDx[i]) / diffQ; + else + A[i][stind] = 0; + } + U0[i][i] = x[i1] - q[i0] * A[i][i]; + U1[i][i] = 2 * x[i2] - q[i1] * A[i][i]; + tx[i] = t; + + if (change[i] == TRUE) { + nTime[i] = t; + } else { + coeff[1] = q[i1] - x[i1]; + coeff[2] = -x[i2]; + if (A[i][i] < 0 && self && x[i2] != 0) { + nTime[i] = t + fabs((q[i1] - x[i1]) / x[i2]); + } else { + nTime[i] = INF; + } + coeff[0] = q[i0] + 2 * lqu[i] - x[i0]; + timeaux = t + minPosRoot(coeff, 2); + if (timeaux < nTime[i]) { + nTime[i] = timeaux; + } + coeff[0] = q[i0] - 2 * lqu[i] - x[i0]; + timeaux = t + minPosRoot(coeff, 2); + if (timeaux < nTime[i]) { + nTime[i] = timeaux; + } + + if (q[i0] * q[i1] < 0 && fabs(q[i0]) > 10 * lqu[i]) { + timeaux = -q[i0] / q[i1] - 2 * fabs(lqu[i] / q[i1]); + if (nTime[i] > t + timeaux) nTime[i] = t + timeaux; + } + double err1 = q[i0] - x[i0] + coeff[1] * (nTime[i] - t) / 2 + coeff[2] * pow((nTime[i] - t) / 2, 2); + if (fabs(err1) > 3 * fabs(lqu[i])) { + nTime[i] = t + quantizer->state->finTime * quantizer->state->minStep; + } + } +} + +#ifdef QSS_PARALLEL +void QSS_TEST_PAR_recomputeNextTimes(QA_quantizer quantizer, int vars, int *inf, double t, double *nTime, double *x, double *lqu, double *q) +#else +void QSS_TEST_recomputeNextTimes(QA_quantizer quantizer, int vars, int *inf, double t, double *nTime, double *x, double *lqu, double *q) +#endif +{ + int i; +#ifdef QSS_PARALLEL + int *map = quantizer->state->qMap; +#endif + for (i = 0; i < vars; i++) { +#ifdef QSS_PARALLEL + if (map[inf[i]] != NOT_ASSIGNED) { +#endif + QSS_TEST_recomputeNextTime(quantizer, inf[i], t, nTime, x, lqu, q); +#ifdef QSS_PARALLEL + } +#endif + } +} + +#ifdef QSS_PARALLEL +void QSS_TEST_PAR_nextTime(QA_quantizer quantizer, int i, double t, double *nTime, double *x, double *lqu) +#else +void QSS_TEST_nextTime(QA_quantizer quantizer, int i, double t, double *nTime, double *x, double *lqu) +#endif +{ + int i2 = i * 3 + 2; + if (x[i2] == 0) + nTime[i] = INF; + else + nTime[i] = t + sqrt(fabs(lqu[i] / x[i2])); +} + +#ifdef QSS_PARALLEL +void QSS_TEST_PAR_solve_single(QA_quantizer quantizer, int i, double *x, double *q, double *lqu, double *h2) +#else +void QSS_TEST_solve_single(QA_quantizer quantizer, int i, double *x, double *q, double *lqu, double *h2) +#endif +{ + int i0 = i * 3, i1 = i0 + 1, i2 = i1 + 1; + double t = quantizer->state->lSimTime->time; + double **U0 = quantizer->state->U0; + double **U1 = quantizer->state->U1; + double **A = quantizer->state->A; + double h, ddx; + + if (A[i][i] < 0) { + ddx = x[i2]; + if (ddx == 0) { + ddx = A[i][i] * A[i][i] * q[i0] + A[i][i] * U0[i][i] + U1[i][i]; + if (ddx == 0) ddx = 1e-40; + } + h = (quantizer->state->finTime - t); + q[i0] = ((x[i0] + h * U0[i][i] + h * h / 2 * U1[i][i]) * (1 - h * A[i][i]) + (h * h / 2 * A[i][i] - h) * (U0[i][i] + h * U1[i][i])) / + (1 - h * A[i][i] + h * h * A[i][i] * A[i][i] / 2); + + if (fabs(q[i0] - x[i0]) > 2 * lqu[i]) { + h = sqrt(fabs(lqu[i] / ddx)); + q[i0] = ((x[i0] + h * U0[i][i] + h * h / 2 * U1[i][i]) * (1 - h * A[i][i]) + (h * h / 2 * A[i][i] - h) * (U0[i][i] + h * U1[i][i])) / + (1 - h * A[i][i] + h * h * A[i][i] * A[i][i] / 2); + } + while (fabs(q[i0] - x[i0]) > 2 * lqu[i]) { + h = h * sqrt(lqu[i] / fabs(q[i0] - x[i0])); + q[i0] = ((x[i0] + h * U0[i][i] + h * h / 2 * U1[i][i]) * (1 - h * A[i][i]) + (h * h / 2 * A[i][i] - h) * (U0[i][i] + h * U1[i][i])) / + (1 - h * A[i][i] + h * h * A[i][i] * A[i][i] / 2); + } + q[i1] = (A[i][i] * q[i0] + U0[i][i] + h * U1[i][i]) / (1 - h * A[i][i]); + } else { + if (x[i2] != 0) { + if (x[i2] > 0) + q[i0] = x[i0] - lqu[i]; + else + q[i0] = x[i0] + lqu[i]; + h = sqrt(2 * (x[i0] - q[i0]) / x[i2]); + h *= 2; + q[i1] = x[i1] + x[i2] * h; + } else { + q[i0] = x[i0]; + q[i1] = x[i1]; + } + } + *h2 = h; +} + +#ifdef QSS_PARALLEL +void QSS_TEST_PAR_AxB(double A[2][2], double B[2][2], double AB[2][2]) +#else +void QSS_TEST_AxB(double A[2][2], double B[2][2], double AB[2][2]) +#endif +{ + AB[0][0] = A[0][0] * B[0][0] + A[0][1] * B[1][0]; + AB[0][1] = A[0][0] * B[0][1] + A[0][1] * B[1][1]; + AB[1][0] = A[1][0] * B[0][0] + A[1][1] * B[1][0]; + AB[1][1] = A[1][0] * B[0][1] + A[1][1] * B[1][1]; +} + +#ifdef QSS_PARALLEL +void QSS_TEST_PAR_invA(double A[2][2], double invA[2][2]) +#else +void QSS_TEST_invA(double A[2][2], double invA[2][2]) +#endif +{ + double detA = A[0][0] * A[1][1] - A[0][1] * A[1][0]; + if (detA == 0) detA = 1e-10; + invA[0][0] = A[1][1] / detA; + invA[0][1] = -A[0][1] / detA; + invA[1][0] = -A[1][0] / detA; + invA[1][1] = A[0][0] / detA; +} + +#ifdef QSS_PARALLEL +void QSS_TEST_PAR_solver2x2_h(QA_quantizer quantizer, double *x, double *q, double *next, int i, int j, double h, double xj0) +#else +void QSS_TEST_solver2x2_h(QA_quantizer quantizer, double *x, double *q, double *next, int i, int j, double h, double xj0) +#endif +{ + int i0 = 3 * i, i1 = i0 + 1; + int j0 = 3 * j, j1 = j0 + 1; + double **U0 = quantizer->state->U0; + double **U1 = quantizer->state->U1; + double **A = quantizer->state->A; + double Ad[2][2] = {{1 - h * A[i][i], -h * A[i][j]}, {-h * A[j][i], 1 - h * A[j][j]}}; + double An[2][2] = {{A[i][i], A[i][j]}, {A[j][i], A[j][j]}}; + double B[2][2], BA[2][2], AB[2][2], ABA[2][2], invC[2][2]; + QSS_TEST_invA(Ad, B); + QSS_TEST_AxB(B, An, BA); + QSS_TEST_AxB(An, B, AB); + QSS_TEST_AxB(An, BA, ABA); + // C = I + h(BA - A - h/2*ABA) + double C[2][2] = {{1 + h * (BA[0][0] - A[i][i] - h / 2 * ABA[0][0]), h * (BA[0][1] - A[i][j] - h / 2 * ABA[0][1])}, + {h * (BA[1][0] - A[j][i] - h / 2 * ABA[1][0]), 1 + h * (BA[1][1] - A[j][j] - h / 2 * ABA[1][1])}}; + // D = I - B + h/2*AB + double D[2][2] = {{1 - B[0][0] + h / 2 * AB[0][0], -B[0][1] + h / 2 * AB[0][1]}, + {-B[1][0] + h / 2 * AB[1][0], 1 - B[1][1] + h / 2 * AB[1][1]}}; + // E = (1/2*(h*AB + I) - B) + double E[2][2] = {{1 / 2 * (h * AB[0][0] + 1) - B[0][0], h / 2 * AB[0][1] - B[0][1]}, + {h / 2 * AB[1][0] - B[1][0], 1 / 2 * (h * AB[1][1] + 1) - B[1][1]}}; + double Du0[2] = {D[0][0] * U0[i][j] + D[0][1] * U0[j][i], D[1][0] * U0[i][j] + D[1][1] * U0[j][i]}; + double Eu1[2] = {E[0][0] * U1[i][j] + E[0][1] * U1[j][i], E[1][0] * U1[i][j] + E[1][1] * U1[j][i]}; + QSS_TEST_invA(C, invC); + q[i0] = invC[0][0] * (x[i0] + h * Du0[0] + h * h * Eu1[0]) + invC[0][1] * (xj0 + h * Du0[1] + h * h * Eu1[1]); + next[j0] = invC[1][0] * (x[i0] + h * Du0[0] + h * h * Eu1[0]) + invC[1][1] * (xj0 + h * Du0[1] + h * h * Eu1[1]); + q[i1] = B[0][0] * (A[i][i] * q[i0] + A[i][j] * next[j0] + U0[i][j] + h * U1[i][j]) + + B[0][1] * (A[j][i] * q[i0] + A[j][j] * next[j0] + U0[j][i] + h * U1[j][i]); + next[j1] = B[1][0] * (A[i][i] * q[i0] + A[i][j] * next[j0] + U0[i][j] + h * U1[i][j]) + + B[1][1] * (A[j][i] * q[i0] + A[j][j] * next[j0] + U0[j][i] + h * U1[j][i]); +} + +#ifdef QSS_PARALLEL +void QSS_TEST_PAR_old_dx(QA_quantizer quantizer, int i, double t, int nSD, double *x, double *tx) +#else +void QSS_TEST_old_dx(QA_quantizer quantizer, int i, double t, int nSD, double *x, double *tx) +#endif +{ + int m, j, j0, j1, j2; + if (t > 0) { + for (m = 0; m < nSD; m++) { + j = quantizer->state->SD[i][m]; + j0 = j * 3; + j1 = j0 + 1; + j2 = j1 + 1; + quantizer->state->oldDx[j] = x[j1] + (t - tx[j]) * x[j2] * 2; + // tx[j] = t; + } + } +} + +#ifdef QSS_PARALLEL +int QSS_TEST_PAR_find_max_inf(QA_quantizer quantizer, double *q, double *x, int i, double qi1_old) +#else +int QSS_TEST_find_max_inf(QA_quantizer quantizer, double *q, double *x, int i, double qi1_old) +#endif +{ + double t = quantizer->state->lSimTime->time; + double **U0 = quantizer->state->U0; + double **U1 = quantizer->state->U1; + double **A = quantizer->state->A; + double *tx = quantizer->state->tx; + int nSD = quantizer->state->nSD[i]; + int m, j, i1 = 3 * i + 1, j1, j2, jmax = -1; + double elapsed, ddxj, max = 0; + + if (t > 0) { + for (m = 0; m < nSD; m++) { + j = quantizer->state->SD[i][m]; + j1 = 3 * j + 1; + j2 = j1 + 1; + if (j != i && A[i][j] * A[j][i] != 0) { + elapsed = t - tx[j]; + // u0 = U0[j][j] + elapsed * U1[j][j]; + // u1 = U1[j][j] - A[j][i] * qi1_old; + ddxj = A[j][i] * q[i1] + A[j][j] * q[j1] + U1[j][i]; + if (fabs(2 * x[j2] - ddxj) > fabs(2 * x[j2] + ddxj) / 10) { + if (fabs(2 * x[j2] - ddxj) / (fabs(2 * x[j2] + ddxj)) > max) { + max = fabs(2 * x[j2] - ddxj) / (fabs(2 * x[j2] + ddxj)); + jmax = j; + } + } + } + } + } + return jmax; +} + +#ifdef QSS_PARALLEL +void QSS_TEST_PAR_updateQuantizedState(QA_quantizer quantizer, int i, double *q, double *x, double *lqu) +#else +void QSS_TEST_updateQuantizedState(QA_quantizer quantizer, int i, double *q, double *x, double *lqu) +#endif +{ + int i0 = i * 3, i1 = i0 + 1, i2 = i1 + 1; + bool *change = quantizer->state->change; + double *next = quantizer->state->next; + double t = quantizer->state->lSimTime->time; + double **U0 = quantizer->state->U0; + double **U1 = quantizer->state->U1; + double **A = quantizer->state->A; + double *tx = quantizer->state->tx; + int nSD = quantizer->state->nSD[i]; + double elapsed; + + double *c = quantizer->state->cont; + double h; + + elapsed = t - quantizer->state->lSimTime->tq[i]; + quantizer->state->qAux[i] = q[i0] + elapsed * q[i1]; + quantizer->state->oldDx[i] = x[i1]; + elapsed = t - tx[i]; + U0[i][i] = U0[i][i] + elapsed * U1[i][i]; + tx[i] = t; + + if (change[i] == TRUE) { + q[i0] = next[i0]; + q[i1] = next[i1]; + change[i] = FALSE; + } else { + double qi_proy, qj_proy, ddxi, qiold[2], qjold[2]; + int j, j0, j1, j2; + + qiold[0] = q[i0]; + qiold[1] = q[i1]; + QSS_TEST_solve_single(quantizer, i, x, q, lqu, &h); +#ifdef QSS_PARALLEL + j = QSS_TEST_PAR_find_max_inf(quantizer, q, x, i, qiold[1]); +#else + j = QSS_TEST_find_max_inf(quantizer, q, x, i, qiold[1]); +#endif + if (t > 0 && j != -1) { + j0 = j * 3; + j1 = j0 + 1; + j2 = j1 + 1; + if (j != i && A[i][j] * A[j][i] != 0) { + double ddxj; + elapsed = t - quantizer->state->lSimTime->tq[j]; + qj_proy = q[j0] + elapsed * q[j1]; + elapsed = t - tx[j]; + U0[j][j] = U0[j][j] + elapsed * U1[j][j]; + U1[j][i] = U1[j][j] - A[j][i] * qiold[1]; + ddxj = A[j][i] * q[i1] + A[j][j] * q[j1] + U1[j][i]; + if (fabs(2 * x[j2] - ddxj) > fabs(2 * x[j2] + ddxj) / 2 || c[0] == j) { + double xj0 = x[j0] + (t - tx[j]) * x[j1] + (t - tx[j]) * (t - tx[j]) * x[j2]; + qi_proy = quantizer->state->qAux[i]; + qjold[0] = q[j0]; + qjold[1] = q[j1]; + U0[j][i] = U0[j][j] - A[j][i] * qi_proy; + U1[i][j] = U1[i][i] - A[i][j] * q[j1]; + if (ddxj > 0) + q[j0] = xj0 - 2 * lqu[j]; + else + q[j0] = xj0 + 2 * lqu[j]; + q[j1] = (A[j][i] * (q[i0] + h * q[i1]) + A[j][j] * q[j0] + U0[j][i] + h * U1[j][i]) / (1 - h * A[j][j]); + ddxi = A[i][i] * q[i1] + A[i][j] * q[j1] + U1[i][j]; + q[j0] = qjold[0]; + q[j1] = qjold[1]; + double a1 = A[i][i] * A[j][j]; + double a2 = A[i][j] * A[j][i]; + bool dA = FALSE; + double tol = 1e-1; + if (a2 != 0 && a1 / a2 < 1 + tol && a1 / a2 > 1 - tol) dA = TRUE; + if (x[i2] * ddxi <= 0 || fabs(x[i2] * ddxi) < 1e-30 || dA) { + U0[i][j] = U0[i][i] - A[i][j] * qj_proy; + qiold[0] = q[i0]; + qiold[1] = q[i1]; + qjold[0] = q[j0]; + qjold[1] = q[j1]; + next[j0] = 1e20; + next[j1] = 1e20; + + h = (quantizer->state->finTime - t) * 1; + QSS_TEST_solver2x2_h(quantizer, x, q, next, i, j, h, xj0); + if (fabs(q[i0] - x[i0]) < 2 * lqu[i] && fabs(next[j0] - xj0) < 2 * lqu[j]) + change[j] = TRUE; + else { + if (fabs(q[i0] - x[i0]) / lqu[i] > fabs(next[j0] - xj0) / lqu[j]) + h = h / sqrt(2 * fabs(q[i0] - x[i0]) / lqu[i]); + else + h = h / sqrt(2 * fabs(next[j0] - xj0) / lqu[j]); + h /= 2; + QSS_TEST_solver2x2_h(quantizer, x, q, next, i, j, h, xj0); + if (fabs(q[i0] - x[i0]) < 2 * lqu[i] && fabs(next[j0] - xj0) < 2 * lqu[j]) + change[j] = TRUE; + else { + if (ddxi != 0 && ddxj != 0) { + if (fabs(lqu[i] / ddxi) < fabs(lqu[j] / ddxj)) + h = sqrt(2 * fabs(lqu[i] / ddxi)); + else + h = sqrt(2 * fabs(lqu[j] / ddxj)); + } else if (ddxi != 0) + h = sqrt(2 * fabs(lqu[i] / ddxi)); + else if (ddxj != 0) + h = sqrt(2 * fabs(lqu[j] / ddxj)); + h /= 2; + QSS_TEST_solver2x2_h(quantizer, x, q, next, i, j, h, xj0); + if (fabs(q[i0] - x[i0]) < 2 * lqu[i] && fabs(next[j0] - xj0) < 2 * lqu[j]) + change[j] = TRUE; + else { + h /= 10; + QSS_TEST_solver2x2_h(quantizer, x, q, next, i, j, h, xj0); + if (fabs(q[i0] - x[i0]) < 2 * lqu[i] && fabs(next[j0] - xj0) < 2 * lqu[j]) + change[j] = TRUE; + else { + q[i0] = qiold[0]; + q[i1] = qiold[1]; + q[j0] = qjold[0]; + q[j1] = qjold[1]; + next[j0] = 1e20; + next[j1] = 1e20; + } + } + } + } + } + } + } + } + } + QSS_TEST_old_dx(quantizer, i, t, nSD, x, tx); + c[0] = i; +} diff --git a/src/engine/qss/methods/qss_test.h b/src/engine/qss/methods/qss_test.h new file mode 100644 index 00000000..3c5c51d9 --- /dev/null +++ b/src/engine/qss/methods/qss_test.h @@ -0,0 +1,76 @@ +/***************************************************************************** + + This file is part of QSS Solver. + + QSS Solver is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + QSS Solver is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with QSS Solver. If not, see . + + ******************************************************************************/ + +#ifndef QSS_TEST_H_ +#define QSS_TEST_H_ + +#include +#include + +extern void QSS_TEST_recomputeNextTimes(QA_quantizer quantizer, int vars, int *inf, double t, double *nTime, double *x, double *lqu, + double *q) __attribute__((hot)); + +extern void QSS_TEST_recomputeNextTime(QA_quantizer quantizer, int var, double t, double *nTime, double *x, double *lqu, double *q) + __attribute__((hot)); + +extern void QSS_TEST_nextTime(QA_quantizer quantizer, int var, double t, double *nTime, double *x, double *lqu) __attribute__((hot)); + +extern void QSS_TEST_updateQuantizedState(QA_quantizer quantizer, int i, double *q, double *x, double *lqu) __attribute__((hot)); + +extern void QSS_TEST_init(QA_quantizer quantizer, QSS_data simData, QSS_time simTime); + +extern void QSS_TEST_PAR_recomputeNextTimes(QA_quantizer quantizer, int vars, int *inf, double t, double *nTime, double *x, double *lqu, + double *q) __attribute__((hot)); + +extern void QSS_TEST_PAR_recomputeNextTime(QA_quantizer quantizer, int var, double t, double *nTime, double *x, double *lqu, double *q) + __attribute__((hot)); + +extern void QSS_TEST_PAR_nextTime(QA_quantizer quantizer, int var, double t, double *nTime, double *x, double *lqu) __attribute__((hot)); + +double qsolver_liqss2(QA_quantizer quantizer, double *x, double *u, double A, double dQ, double *q); + +void solve_2x2(QA_quantizer quantizer, int i, int j, double xi0, double xj0, double qi0, double *qj, double *h); + +void old_dx(QA_quantizer quantizer, int i, double t, int nSD, double *x, double *tx); + +extern void QSS_TEST_PAR_updateQuantizedState(QA_quantizer quantizer, int i, double *q, double *x, double *lqu) __attribute__((hot)); + +extern void QSS_TEST_PAR_init(QA_quantizer quantizer, QSS_data simData, QSS_time simTime); + +extern void QSS_TEST_PAR_AxB(double A[2][2], double B[2][2], double AB[2][2]); + +extern void QSS_TEST_AxB(double A[2][2], double B[2][2], double AB[2][2]); + +extern void QSS_TEST_PAR_invA(double A[2][2], double invA[2][2]); + +extern void QSS_TEST_invA(double A[2][2], double invA[2][2]); + +extern void QSS_TEST_PAR_solver2x2_h(QA_quantizer quantizer, double *x, double *q, double *next, int i, int j, double h, double xj0); + +extern void QSS_TEST_solver2x2_h(QA_quantizer quantizer, double *x, double *q, double *next, int i, int j, double h, double xj0); + +extern void QSS_TEST_PAR_old_dx(QA_quantizer quantizer, int i, double t, int nSD, double *x, double *tx); + +extern void QSS_TEST_old_dx(QA_quantizer quantizer, int i, double t, int nSD, double *x, double *tx); + +extern void QSS_TEST_PAR_solve_single(QA_quantizer quantizer, int i, double *x, double *q, double *lqu, double *h2); + +extern void QSS_TEST_solve_single(QA_quantizer quantizer, int i, double *x, double *q, double *lqu, double *h2); + +#endif /* QSS_TEST_H_ */ diff --git a/src/engine/qss/parallel/methods/Makefile.include b/src/engine/qss/parallel/methods/Makefile.include index d91504cd..fc48c619 100644 --- a/src/engine/qss/parallel/methods/Makefile.include +++ b/src/engine/qss/parallel/methods/Makefile.include @@ -11,7 +11,8 @@ PAR_METHODS_SRC = $(QSS_METHODS_DIR)/qss.c \ $(QSS_METHODS_DIR)/liqss.c \ $(QSS_METHODS_DIR)/liqss2.c \ $(QSS_METHODS_DIR)/liqss_bdf.c \ - $(QSS_METHODS_DIR)/liqss3.c + $(QSS_METHODS_DIR)/liqss3.c \ + $(QSS_METHODS_DIR)/qss_test.c \ # Objects diff --git a/src/engine/qss/qss_quantizer.c b/src/engine/qss/qss_quantizer.c index 31e85076..a298f7d0 100644 --- a/src/engine/qss/qss_quantizer.c +++ b/src/engine/qss/qss_quantizer.c @@ -32,6 +32,7 @@ #include #include #include +#include QA_quantizerOps QA_QuantizerOps() { @@ -72,6 +73,16 @@ QA_quantizerState QA_QuantizerState() p->finTime = 0; p->lSimTime = NULL; p->qMap = NULL; + p->nSD = NULL; + p->SD = NULL; + p->change = NULL; + p->next = NULL; + p->nTime = NULL; + p->A = NULL; + p->U0 = NULL; + p->U1 = NULL; + p->cont = NULL; + return p; } @@ -144,6 +155,13 @@ QA_quantizer QA_Quantizer(QSS_data simData, QSS_time simTime) QSS4_init(p, simData, simTime); } break; + case SD_QSS_TEST: + if (simData->params->lps > 0) { + QSS_TEST_PAR_init(p, simData, simTime); + } else { + QSS_TEST_init(p, simData, simTime); + } + break; default: return NULL; } diff --git a/src/engine/qss/qss_quantizer.h b/src/engine/qss/qss_quantizer.h index fc22c537..61b636b5 100644 --- a/src/engine/qss/qss_quantizer.h +++ b/src/engine/qss/qss_quantizer.h @@ -23,176 +23,84 @@ #include #include -/** - * - */ typedef struct QA_quantizerOps_ *QA_quantizerOps; -/** - * - */ typedef struct QA_quantizerState_ *QA_quantizerState; -/** - * - */ typedef struct QA_quantizer_ *QA_quantizer; -/** - * - * @param - * @param - * @param - * @param - * @param - * @param - * @param - */ typedef void (*QA_recomputeNextTimesFn)(QA_quantizer, int, int *, double, double *, double *, double *, double *); -/** - * - * @param - * @param - * @param - * @param - * @param - * @param - */ typedef void (*QA_recomputeNextTimeFn)(QA_quantizer, int, double, double *, double *, double *, double *); -/** - * - * @param - * @param - * @param - * @param - * @param - */ typedef void (*QA_nextTimeFn)(QA_quantizer, int, double, double *, double *, double *); -/** - * - * @param - * @param - * @param - * @param - */ typedef void (*QA_updateQuantizedStateFn)(QA_quantizer, int, double *, double *, double *); -/** - * - */ + struct QA_quantizerState_ { int order; int xOrder; - int *flag2; //!< - int *flag3; //!< - int *flag4; //!< + int *flag2; + int *flag3; + int *flag4; int *BDF; int *nSZ; int **SZ; - double *dq; //!< - double *a; //!< - double *oldDx; //!< - double *qAux; //!< - double *u0; //!< - double *u1; //!< - double *u2; //!< + double *dq; + double *a; + double *oldDx; + double *qAux; + double *u0; + double *u1; + double *u2; double *tx; - double *lt; //!< - double *ltq; //!< - double *lquOld; //!< - double *qinf; //!< - double *qsup; //!< - double *simTime; //!< - double minStep; //!< - double finTime; //!< - QSS_time lSimTime; //!< - QSS_idxMap qMap; //!< + double *lt; + double *ltq; + double *lquOld; + double *qinf; + double *qsup; + double *simTime; + double minStep; + double finTime; + QSS_time lSimTime; + QSS_idxMap qMap; + + // QSS TEST method added state. + int *nSD; + int **SD; + bool *change; + double *next; + double *nTime; + double **A; + double **U0; + double **U1; + double *cont; }; -/** - * - */ + struct QA_quantizerOps_ { - QA_recomputeNextTimesFn recomputeNextTimes; //!< - QA_recomputeNextTimeFn recomputeNextTime; //!< - QA_nextTimeFn nextTime; //!< - QA_updateQuantizedStateFn updateQuantizedState; //!< + QA_recomputeNextTimesFn recomputeNextTimes; + QA_recomputeNextTimeFn recomputeNextTime; + QA_nextTimeFn nextTime; + QA_updateQuantizedStateFn updateQuantizedState; }; -/** - * - */ + struct QA_quantizer_ { - QA_quantizerState state; //!< - QA_quantizerOps ops; //!< + QA_quantizerState state; + QA_quantizerOps ops; }; -/** - * - * @return - */ + QA_quantizerOps QA_QuantizerOps(); -/** - * - * @return - */ + QA_quantizerState QA_QuantizerState(); -/** - * - * @param simData - * @param simTime - * @return - */ + QA_quantizer QA_Quantizer(QSS_data simData, QSS_time simTime); -/** - * - * @param state - */ + void QA_freeQuantizerState(QA_quantizerState state); -/** - * - * @param ops - */ + void QA_freeQuantizerOps(QA_quantizerOps ops); -/** - * - * @param q - */ + void QA_freeQuantizer(QA_quantizer quantizer); -/** - * - * @param quantizer - * @param vars - * @param inf - * @param t - * @param nTime - * @param x - * @param lqu - * @param q - */ + void QA_recomputeNextTimes(QA_quantizer quantizer, int vars, int *inf, double t, double *nTime, double *x, double *lqu, double *q); -/** - * - * @param quantizer - * @param var - * @param t - * @param nTime - * @param x - * @param lqu - * @param q - */ + void QA_recomputeNextTime(QA_quantizer quantizer, int var, double t, double *nTime, double *x, double *lqu, double *q); -/** - * - * @param quantizer - * @param var - * @param t - * @param nTime - * @param x - * @param lqu - */ + void QA_nextTime(QA_quantizer quantizer, int var, double t, double *nTime, double *x, double *lqu); -/** - * - * @param quantizer - * @param x - * @param q - * @param x - * @param lqu - */ + void QA_updateQuantizedState(QA_quantizer quantizer, int var, double *q, double *x, double *lqu); #endif /* QSS_QUANTIZER_H_ */