diff --git a/lib/algo/IncrementalDMD.cpp b/lib/algo/IncrementalDMD.cpp index 32ab8c84b..9b3836e17 100644 --- a/lib/algo/IncrementalDMD.cpp +++ b/lib/algo/IncrementalDMD.cpp @@ -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; } @@ -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) { diff --git a/lib/linalg/svd/IncrementalSVDBrand.cpp b/lib/linalg/svd/IncrementalSVDBrand.cpp index f5463f40b..0a987b1fc 100644 --- a/lib/linalg/svd/IncrementalSVDBrand.cpp +++ b/lib/linalg/svd/IncrementalSVDBrand.cpp @@ -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) {