14 #include "RandomizedSVD.h"
17 #include "linalg/scalapack_wrapper.h"
18 #include "utils/mpi_utils.h"
29 RandomizedSVD::RandomizedSVD(
32 d_subspace_dim(options.randomized_subspace_dim) {
33 srand(options.random_seed);
37 RandomizedSVD::computeSVD()
39 d_samples->n = d_num_samples;
42 const int num_rows = d_total_dim;
43 const int num_cols = d_num_samples;
44 if (d_subspace_dim < 1 || d_subspace_dim > std::min(num_rows, num_cols)) {
45 d_subspace_dim = std::min(num_rows, num_cols);
50 Matrix* snapshot_matrix;
51 if (num_rows > num_cols) {
52 snapshot_matrix =
new Matrix(d_dims[
static_cast<unsigned>(d_rank)],
54 for (
int rank = 0; rank < d_num_procs; ++rank) {
57 gather_transposed_block(&snapshot_matrix->item(0, 0), d_samples.get(),
58 d_istarts[
static_cast<unsigned>(rank)] + 1, 1,
59 d_dims[
static_cast<unsigned>(rank)], num_cols,
64 std::vector<int> snapshot_transpose_row_offset;
65 const int num_transposed_rows =
split_dimension(num_cols, MPI_COMM_WORLD);
67 snapshot_transpose_row_offset,
69 CAROM_VERIFY(num_cols == num_cols_check);
71 snapshot_matrix =
new Matrix(num_transposed_rows,
74 for (
int rank = 0; rank < d_num_procs; ++rank) {
75 gather_block(&snapshot_matrix->item(0, 0), d_samples.get(),
76 1, snapshot_transpose_row_offset[rank] + 1,
77 num_rows, snapshot_transpose_row_offset[rank + 1] -
78 snapshot_transpose_row_offset[rank],
83 int snapshot_matrix_distributed_rows = std::max(num_rows, num_cols);
89 if (d_debug_algorithm) {
90 rand_mat =
new Matrix(snapshot_matrix->numColumns(), d_subspace_dim,
false);
91 for (
int i = 0; i < snapshot_matrix->numColumns(); i++) {
92 for (
int j = 0; j < d_subspace_dim; j++) {
93 rand_mat->item(i, j) = (i == j);
98 rand_mat =
new Matrix(snapshot_matrix->numColumns(), d_subspace_dim,
false,
103 Matrix* rand_proj = snapshot_matrix->mult(rand_mat);
104 int rand_proj_rows = rand_proj->numRows();
107 Matrix* Q = rand_proj->qr_factorize();
110 Matrix* svd_input_mat = Q->transposeMult(snapshot_matrix);
111 int svd_input_mat_distributed_rows = svd_input_mat->numDistributedRows();
113 SLPK_Matrix svd_input;
114 initialize_matrix(&svd_input, svd_input_mat->numColumns(),
115 svd_input_mat->numDistributedRows(),
116 d_npcol, d_nprow, d_blocksize, d_blocksize);
117 scatter_block(&svd_input, 1, 1,
118 svd_input_mat->getData(),
119 svd_input_mat->numColumns(), svd_input_mat->numDistributedRows(),0);
120 delete svd_input_mat;
121 delete snapshot_matrix;
124 svd_init(d_factorizer.get(), &svd_input);
126 d_factorizer->dov = 1;
128 factorize(d_factorizer.get());
129 free_matrix_data(&svd_input);
132 int sigma_cutoff = 0, hard_cutoff = num_cols;
133 if (d_singular_value_tol == 0) {
134 sigma_cutoff = std::numeric_limits<int>::max();
136 for (
int i = 0; i < num_cols; ++i) {
137 if (d_factorizer->S[i] / d_factorizer->S[0] > d_singular_value_tol) {
144 if (d_max_basis_dimension != -1 && d_max_basis_dimension < hard_cutoff) {
145 hard_cutoff = d_max_basis_dimension;
147 int ncolumns = hard_cutoff < sigma_cutoff ? hard_cutoff : sigma_cutoff;
148 CAROM_VERIFY(ncolumns >= 0);
149 ncolumns = std::min(ncolumns, d_subspace_dim);
152 d_S =
new Vector(ncolumns,
false);
154 CAROM_VERIFY(ncolumns >= 0);
155 unsigned nc =
static_cast<unsigned>(ncolumns);
156 memset(&d_S->item(0), 0, nc*
sizeof(
double));
159 d_basis =
new Matrix(svd_input_mat_distributed_rows, ncolumns,
false);
160 d_basis_right =
new Matrix(d_subspace_dim, ncolumns,
false);
162 for (
int rank = 0; rank < d_num_procs; ++rank) {
164 gather_block(&d_basis->item(0, 0), d_factorizer->V, 1, 1,
165 ncolumns, d_subspace_dim, rank);
168 gather_transposed_block(&d_basis_right->item(0, 0), d_factorizer->U,
170 1, svd_input_mat_distributed_rows,
174 for (
int i = 0; i < ncolumns; ++i)
175 d_S->item(i) = d_factorizer->S[
static_cast<unsigned>(i)];
178 Matrix* d_new_basis = Q->mult(d_basis);
180 d_basis = d_new_basis;
182 if (num_rows <= num_cols) {
183 Matrix* temp = d_basis;
184 d_basis = d_basis_right;
185 d_basis_right = temp;
188 if (d_basis->numRows() == d_total_dim)
191 printf(
"WARNING: RandomizedSVD::d_basis has a subspace dimension %d, different from sample dimension %d\n",
192 d_basis->numRows(), num_rows);
195 d_basis->distribute(local_rows);
196 d_basis_right->gather();
199 d_basis_is_current =
true;
201 release_context(&svd_input);
203 if (d_debug_algorithm) {
205 printf(
"Distribution of sampler's A and U:\n");
207 print_debug_info(d_samples.get());
208 MPI_Barrier(MPI_COMM_WORLD);
211 printf(
"Distribution of sampler's V:\n");
213 print_debug_info(d_factorizer->V);
214 MPI_Barrier(MPI_COMM_WORLD);
217 printf(
"Computed singular values: ");
218 for (
int i = 0; i < ncolumns; ++i)
219 printf(
"%8.4E ", d_factorizer->S[i]);
int split_dimension(const int dim, const MPI_Comm &comm)
Distribute the global size dim into MPI processes as equally as possible.
int get_global_offsets(const int local_dim, std::vector< int > &offsets, const MPI_Comm &comm)
Save integer offsets for each MPI rank under MPI communicator comm, where each rank as the size of lo...