Skip to content

Commit

Permalink
stylize
Browse files Browse the repository at this point in the history
  • Loading branch information
swsuh28 committed Apr 9, 2024
1 parent 06b4d6b commit 9442d25
Show file tree
Hide file tree
Showing 2 changed files with 79 additions and 79 deletions.
156 changes: 78 additions & 78 deletions lib/algo/IncrementalDMD.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -134,54 +134,54 @@ IncrementalDMD::updateDMD(const Matrix* f_snapshots)

int num_samples = svd->getNumSamples();

if (num_samples == 1)
{
Matrix* d_basis = new Matrix(*(svd->getSpatialBasis()));
Matrix* d_basis_right = new Matrix(*(svd->getTemporalBasis()));
Vector* d_sv = new Vector(*(svd->getSingularValues()));
Matrix* d_S_inv = new Matrix(1, 1, false);
d_S_inv->item(0, 0) = 1/d_sv->item(0);

Matrix* f_snapshots_out = new Matrix(d_snapshots.back()->getData(),
d_dim, 1, true);
Matrix* br_Sinv = d_basis_right->mult(d_S_inv);
Matrix* f_br_Sinv = f_snapshots_out->mult(br_Sinv);

d_A_tilde = d_basis->transposeMult(f_br_Sinv);

delete br_Sinv;
delete f_br_Sinv;
delete f_snapshots_out;

delete d_basis;
delete d_basis_right;
delete d_sv;
delete d_S_inv;
}
else
{
Vector u_vec(u_in, svd->d_dim, true);
Vector e_proj(u_in, svd->d_dim, true);

Vector* UTe = svd->d_U_pre->transposeMult(e_proj);
Vector* UUTe = svd->d_U_pre->mult(UTe);
e_proj -= *UUTe;
delete UTe;
delete UUTe;
UTe = svd->d_U_pre->transposeMult(e_proj);
UUTe = svd->d_U_pre->mult(UTe);
e_proj -= *UUTe;
delete UTe;
delete UUTe;
double k = e_proj.inner_product(e_proj);
if (k <= 0) k = 0;
else k = sqrt(k);

if ( k > svd->d_linearity_tol ){
if (num_samples == 1)
{
Matrix* d_basis = new Matrix(*(svd->getSpatialBasis()));
Matrix* d_basis_right = new Matrix(*(svd->getTemporalBasis()));
Vector* d_sv = new Vector(*(svd->getSingularValues()));
Matrix* d_S_inv = new Matrix(1, 1, false);
d_S_inv->item(0, 0) = 1/d_sv->item(0);

Matrix* f_snapshots_out = new Matrix(d_snapshots.back()->getData(),
d_dim, 1, true);
Matrix* br_Sinv = d_basis_right->mult(d_S_inv);
Matrix* f_br_Sinv = f_snapshots_out->mult(br_Sinv);

d_A_tilde = d_basis->transposeMult(f_br_Sinv);

delete br_Sinv;
delete f_br_Sinv;
delete f_snapshots_out;

delete d_basis;
delete d_basis_right;
delete d_sv;
delete d_S_inv;
}
else
{
Vector u_vec(u_in, svd->d_dim, true);
Vector e_proj(u_in, svd->d_dim, true);

Vector* UTe = svd->d_U_pre->transposeMult(e_proj);
Vector* UUTe = svd->d_U_pre->mult(UTe);
e_proj -= *UUTe;

delete UTe;
delete UUTe;

UTe = svd->d_U_pre->transposeMult(e_proj);
UUTe = svd->d_U_pre->mult(UTe);
e_proj -= *UUTe;

delete UTe;
delete UUTe;

double k = e_proj.inner_product(e_proj);
if (k <= 0) k = 0;
else k = sqrt(k);

if ( k > svd->d_linearity_tol ) {
if (d_rank == 0) {
std::cout << "Added linearly independent sample" << std::endl;
}
Expand Down Expand Up @@ -237,45 +237,45 @@ IncrementalDMD::updateDMD(const Matrix* f_snapshots)
delete u_new;
delete d_A_tilde_tmp;
}
//}
else
{
if (d_rank == 0) {
std::cout << "Added linearly dependent sample" << std::endl;
}
//}
else
{
if (d_rank == 0) {
std::cout << "Added linearly dependent sample" << std::endl;
}

Matrix* d_A_tilde_tmp = new Matrix(num_samples, num_samples+1, false);
Vector* u_new = new Vector(d_snapshots[num_snapshots-1]->getData(),
d_dim, true);
Matrix* d_A_tilde_tmp = new Matrix(num_samples, num_samples+1, false);
Vector* u_new = new Vector(d_snapshots[num_snapshots-1]->getData(),
d_dim, true);

Vector* UTu = new Vector(num_samples, false);

svd->d_U_pre->transposeMult(*u_new, UTu);
Vector* UpTUTu = svd->d_Up_pre->transposeMult(UTu);
Vector* UTu = new Vector(num_samples, false);

for (int i = 0; i < num_samples; i++) {
for (int j = 0; j < num_samples; j++) {
d_A_tilde_tmp->item(i, j) = d_A_tilde->item(i, j) * svd->d_S_pre->item(j);
svd->d_U_pre->transposeMult(*u_new, UTu);
Vector* UpTUTu = svd->d_Up_pre->transposeMult(UTu);

for (int i = 0; i < num_samples; i++) {
for (int j = 0; j < num_samples; j++) {
d_A_tilde_tmp->item(i, j) = d_A_tilde->item(i, j) * svd->d_S_pre->item(j);
}
d_A_tilde_tmp->item(i, num_samples) = UpTUTu->item(i);
}
d_A_tilde_tmp->item(i, num_samples) = UpTUTu->item(i);
}

Matrix* WSinv = svd->d_Wq->mult(svd->d_Sq_inv);
Matrix* AWSinv = d_A_tilde_tmp->mult(WSinv);
Matrix* WSinv = svd->d_Wq->mult(svd->d_Sq_inv);
Matrix* AWSinv = d_A_tilde_tmp->mult(WSinv);

Matrix* d_A_tilde_new = svd->d_Uq->transposeMult(AWSinv);
Matrix* d_A_tilde_new = svd->d_Uq->transposeMult(AWSinv);

delete UTu;
delete UpTUTu;
delete WSinv;
delete AWSinv;
delete UTu;
delete UpTUTu;
delete WSinv;
delete AWSinv;

delete d_A_tilde;
d_A_tilde = d_A_tilde_new;
delete d_A_tilde;
d_A_tilde = d_A_tilde_new;

delete u_new;
delete d_A_tilde_tmp;
}
delete u_new;
delete d_A_tilde_tmp;
}
}

if (d_rank == 0) {
Expand Down
2 changes: 1 addition & 1 deletion lib/linalg/svd/IncrementalSVDBrand.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -579,7 +579,7 @@ IncrementalSVDBrand::addNewSample(
if (d_update_right_SV) {
new_d_Wp = new Matrix(d_num_samples+1, d_num_samples+1, false);
new_d_Wp_inv = new Matrix(d_num_samples+1, d_num_samples+1, false);

d_W->item(d_num_rows_of_W, d_num_samples) = 1.0;

for (int row = 0; row < d_num_samples; ++row) {
Expand Down

0 comments on commit 9442d25

Please sign in to comment.