<commit_before>#include "../../lower_triangular_solve.hxx"
#include "../../SDP_Solver.hxx"
#include "../../../Timers.hxx"

// Compute the quantities needed to solve the Schur complement
// equation
//
// {{S, -B}, {B^T, 0}} . {dx, dy} = {r, s}
//
// (where S = SchurComplement, B = FreeVarMatrix), using the method
// described in the manual:
//
// - Compute S using BilinearPairingsXInv and BilinearPairingsY.
//
// - Compute the Cholesky decomposition S' = L' L'^T.
//
// - Form B' = (B U) and compute
//
//   - SchurOffDiagonal = L'^{-1} B
//   - L'^{-1} U
//   - Q = (L'^{-1} B')^T (L'^{-1} B') - {{0, 0}, {0, 1}}
//
// - Compute the LU decomposition of Q.
//
// This data is sufficient to efficiently solve the above equation for
// a given r,s.
//
// Inputs:
// - BilinearPairingsXInv, BilinearPairingsY (these are members of
//   SDPSolver, but we include them as arguments to emphasize that
//   they must be computed first)
// Workspace (members of SDPSolver which are modified by this method
// and not used later):
// - SchurComplement
// Outputs (members of SDPSolver which are modified by this method and
// used later):
// - SchurComplementCholesky
// - SchurOffDiagonal
// - Q, Qpivots
//

void compute_schur_complement(
  const SDP &sdp, const Block_Diagonal_Matrix &bilinear_pairings_X_inv,
  const Block_Diagonal_Matrix &bilinear_pairings_Y,
  Block_Diagonal_Matrix &schur_complement);

void SDP_Solver::initialize_schur_complement_solver(
  const Block_Diagonal_Matrix &bilinear_pairings_X_inv,
  const Block_Diagonal_Matrix &bilinear_pairings_Y,
  const std::vector<size_t> &block_dims,
  Block_Diagonal_Matrix &schur_complement_cholesky, Matrix &schur_off_diagonal,
  Block_Matrix &schur_off_diagonal_elemental, Matrix &Q,
  El::DistMatrix<El::BigFloat> &Q_elemental)
{
  // The Schur complement matrix S: a Block_Diagonal_Matrix with one
  // block for each 0 <= j < J.  SchurComplement.blocks[j] has dimension
  // (d_j+1)*m_j*(m_j+1)/2
  //
  Block_Diagonal_Matrix schur_complement(block_dims);

  compute_schur_complement(sdp, bilinear_pairings_X_inv, bilinear_pairings_Y,
                           schur_complement);

  // compute SchurComplementCholesky = L', where
  //
  //   L' L'^T = S'
  //
  timers["initializeSchurComplementSolver.choleskyDecomposition"].resume();
  cholesky_decomposition(schur_complement, schur_complement_cholesky);
  timers["initializeSchurComplementSolver.choleskyDecomposition"].stop();

  // SchurOffDiagonal = L'^{-1} FreeVarMatrix
  schur_off_diagonal = sdp.free_var_matrix;
  schur_off_diagonal_elemental = sdp.free_var_matrix_elemental;
  timers["initializeSchurComplementSolver.blockMatrixLowerTriangularSolve"]
    .resume();
  lower_triangular_solve(schur_complement_cholesky, schur_off_diagonal);
  lower_triangular_solve(schur_complement_cholesky,
                         schur_off_diagonal_elemental);
  timers["initializeSchurComplementSolver.blockMatrixLowerTriangularSolve"]
    .stop();

  // total number of columns in the off-diagonal part B' = (B U)
  // (currently just B; will accumulate the rest shortly)
  int off_diagonal_columns = schur_off_diagonal.cols;

  // Next, we compute
  //
  //   Q = (L'^{-1} B')^T (L'^{-1} B') - {{0, 0}, {0, 1}}
  //
  // Where B' = (B U).  We think of Q as containing four blocks called
  // Upper/Lower-Left/Right.

  timers["initializeSchurComplementSolver.Qcomputation"].resume();
  // Set the dimensions of Q
  Q.resize(off_diagonal_columns, off_diagonal_columns);
  Q.set_zero();

  // Here, SchurOffDiagonal = L'^{-1} B.
  //
  // UpperLeft(Q) = SchurOffDiagonal^T SchurOffDiagonal
  matrix_square_into_block(schur_off_diagonal, Q, 0, 0);

  // const El::Grid grid(El::mpi::COMM_WORLD);
  // El::DistMatrix<El::BigFloat> schur_off_diagonal_elemental(grid),
  //   Q_elemental(grid);
  // schur_off_diagonal_elemental.Resize(schur_off_diagonal.rows,
  //                                     schur_off_diagonal.cols);
  // Zeros(Q_elemental, schur_off_diagonal.cols, schur_off_diagonal.cols);

  // const El::Int local_height(schur_off_diagonal_elemental.LocalHeight()),
  //   local_width(schur_off_diagonal_elemental.LocalWidth());

  // for(size_t row = 0; row < local_height; ++row)
  //   {
  //     El::Int global_row(schur_off_diagonal_elemental.GlobalRow(row));
  //     for(size_t column = 0; column < local_width; ++column)
  //       {
  //         El::Int global_column(
  //           schur_off_diagonal_elemental.GlobalCol(column));

  //         std::stringstream ss;
  //         ss.precision(1024 * 0.31 + 5);
  //         ss << schur_off_diagonal.elt(global_row, global_column);
  //         schur_off_diagonal_elemental.SetLocal(row, column,
  //                                               El::BigFloat(ss.str(), 10));
  //       }
  //   }
  // El::Syrk(El::UPPER, El::TRANSPOSE, El::BigFloat(1),
  //          schur_off_diagonal_elemental, El::BigFloat(1), Q_elemental);

  // const El::Int Q_height(Q_elemental.LocalHeight()),
  //   Q_width(Q_elemental.LocalWidth());

  // for(size_t row = 0; row < Q_height; ++row)
  //   {
  //     El::Int global_row(Q_elemental.GlobalRow(row));
  //     for(size_t column = 0; column < Q_width; ++column)
  //       {
  //         El::Int global_column(Q_elemental.GlobalCol(column));

  //         std::stringstream ss;
  //         ss.precision(1024 * 0.31 + 5);
  //         ss << Q_elemental.GetLocal(row, column);
  //         Q.elt(global_row, global_column) = ss.str();
  //       }
  //   }

  // std::cout.precision(1024 * 0.31 + 5);
  // std::cout << "Qe: "
  //           // << local_height << " " << local_width << " "
  //           << Q_height << " " << Q_width << " " << Q.rows << " " << Q.cols
  //           << " " << Q_elemental.GetLocal(0, 0) << " "
  //           << "\n";
  // std::cout << "Qe: "
  //           // << local_height << " " << local_width << " "
  //           << Q_height << " " << Q_width << " " << Q.rows << " " << Q.cols
  //           << " " << Q.elt(0, 0) << " "
  //           << "\n";
  // std::cout << El::mpfr::Precision() << " "
  //           << schur_off_diagonal.elt(0,0).get_prec() << " "
  //           << Q.elt(0, 0).get_prec() << " "
  //            << "\n";

  // // std::cout << "Qm: " << Q_height << " " << Q_width << " " << Q.rows << "
  // "
  // //           << Q.cols << " " << Q.elt(0, 0) << " "
  // //           << "\n";

  // UpperRight(Q) = LowerLeft(Q)^T
  for(size_t c = 0; c < schur_off_diagonal.cols; c++)
    for(size_t r = schur_off_diagonal.cols; r < Q.rows; r++)
      {
        Q.elt(c, r) = Q.elt(r, c);
      }
  timers["initializeSchurComplementSolver.Qcomputation"].stop();

  // Resize Qpivots appropriately and compute the LU decomposition of Q

  timers["initializeSchurComplementSolver.LUDecomposition"].resume();
  timers["LUDecomposition.actualLU"].resume();

  Matrix Q_temp(Q);
  cholesky_decomposition(Q_temp, Q);
  timers["LUDecomposition.actualLU"].stop();
  timers["initializeSchurComplementSolver.LUDecomposition"].stop();
}<commit_msg>Compute the elemental version of schur_off_diagonal^2<commit_after>#include "../../lower_triangular_solve.hxx"
#include "../../SDP_Solver.hxx"
#include "../../../Timers.hxx"

// Compute the quantities needed to solve the Schur complement
// equation
//
// {{S, -B}, {B^T, 0}} . {dx, dy} = {r, s}
//
// (where S = SchurComplement, B = FreeVarMatrix), using the method
// described in the manual:
//
// - Compute S using BilinearPairingsXInv and BilinearPairingsY.
//
// - Compute the Cholesky decomposition S' = L' L'^T.
//
// - Form B' = (B U) and compute
//
//   - SchurOffDiagonal = L'^{-1} B
//   - L'^{-1} U
//   - Q = (L'^{-1} B')^T (L'^{-1} B') - {{0, 0}, {0, 1}}
//
// - Compute the LU decomposition of Q.
//
// This data is sufficient to efficiently solve the above equation for
// a given r,s.
//
// Inputs:
// - BilinearPairingsXInv, BilinearPairingsY (these are members of
//   SDPSolver, but we include them as arguments to emphasize that
//   they must be computed first)
// Workspace (members of SDPSolver which are modified by this method
// and not used later):
// - SchurComplement
// Outputs (members of SDPSolver which are modified by this method and
// used later):
// - SchurComplementCholesky
// - SchurOffDiagonal
//

void compute_schur_complement(
  const SDP &sdp, const Block_Diagonal_Matrix &bilinear_pairings_X_inv,
  const Block_Diagonal_Matrix &bilinear_pairings_Y,
  Block_Diagonal_Matrix &schur_complement);

void SDP_Solver::initialize_schur_complement_solver(
  const Block_Diagonal_Matrix &bilinear_pairings_X_inv,
  const Block_Diagonal_Matrix &bilinear_pairings_Y,
  const std::vector<size_t> &block_dims,
  Block_Diagonal_Matrix &schur_complement_cholesky, Matrix &schur_off_diagonal,
  Block_Matrix &schur_off_diagonal_block, Matrix &Q,
  El::DistMatrix<El::BigFloat> &Q_elemental)
{
  // The Schur complement matrix S: a Block_Diagonal_Matrix with one
  // block for each 0 <= j < J.  SchurComplement.blocks[j] has dimension
  // (d_j+1)*m_j*(m_j+1)/2
  //
  Block_Diagonal_Matrix schur_complement(block_dims);

  compute_schur_complement(sdp, bilinear_pairings_X_inv, bilinear_pairings_Y,
                           schur_complement);

  // compute SchurComplementCholesky = L', where
  //
  //   L' L'^T = S'
  //
  timers["initializeSchurComplementSolver.choleskyDecomposition"].resume();
  cholesky_decomposition(schur_complement, schur_complement_cholesky);
  timers["initializeSchurComplementSolver.choleskyDecomposition"].stop();

  // SchurOffDiagonal = L'^{-1} FreeVarMatrix
  schur_off_diagonal = sdp.free_var_matrix;
  schur_off_diagonal_block = sdp.free_var_matrix_elemental;
  timers["initializeSchurComplementSolver.blockMatrixLowerTriangularSolve"]
    .resume();
  lower_triangular_solve(schur_complement_cholesky, schur_off_diagonal);
  lower_triangular_solve(schur_complement_cholesky, schur_off_diagonal_block);
  timers["initializeSchurComplementSolver.blockMatrixLowerTriangularSolve"]
    .stop();

  // total number of columns in the off-diagonal part B' = (B U)
  // (currently just B; will accumulate the rest shortly)
  int off_diagonal_columns = schur_off_diagonal.cols;

  // Next, we compute
  //
  //   Q = (L'^{-1} B')^T (L'^{-1} B') - {{0, 0}, {0, 1}}
  //
  // Where B' = (B U).  We think of Q as containing four blocks called
  // Upper/Lower-Left/Right.

  timers["initializeSchurComplementSolver.Qcomputation"].resume();
  // Set the dimensions of Q
  Q.resize(off_diagonal_columns, off_diagonal_columns);
  Q.set_zero();

  const size_t Q_size(schur_off_diagonal_block.blocks.at(0).Width());
  Zeros(Q_elemental, Q_size, Q_size);

  size_t schur_off_diagonal_total_rows(0);
  for(auto &block : schur_off_diagonal_block.blocks)
    {
      schur_off_diagonal_total_rows += block.Height();
    }
  El::DistMatrix<El::BigFloat> schur_off_diagonal_dist(
    schur_off_diagonal_total_rows, Q_size);

  size_t row_offset(0);
  for(auto &block : schur_off_diagonal_block.blocks)
    {
      El::DistMatrix<El::BigFloat> sub_matrix(
        El::View(schur_off_diagonal_dist, row_offset, 0, block.Height(),
                 block.Width()));
      El::Copy(block, sub_matrix);
      row_offset += block.Height();
    }

  // Here, SchurOffDiagonal = L'^{-1} B.
  //
  // UpperLeft(Q) = SchurOffDiagonal^T SchurOffDiagonal
  matrix_square_into_block(schur_off_diagonal, Q, 0, 0);

  // UpperRight(Q) = LowerLeft(Q)^T
  for(size_t c = 0; c < schur_off_diagonal.cols; c++)
    for(size_t r = schur_off_diagonal.cols; r < Q.rows; r++)
      {
        Q.elt(c, r) = Q.elt(r, c);
      }
  
  El::Syrk(El::UpperOrLowerNS::UPPER, El::OrientationNS::TRANSPOSE,
           El::BigFloat(1), schur_off_diagonal_dist, El::BigFloat(1),
           Q_elemental);
  El::MakeSymmetric(El::UpperOrLower::UPPER, Q_elemental);
  timers["initializeSchurComplementSolver.Qcomputation"].stop();

  timers["initializeSchurComplementSolver.LUDecomposition"].resume();
  timers["LUDecomposition.actualLU"].resume();

  Matrix Q_temp(Q);
  cholesky_decomposition(Q_temp, Q);
  timers["LUDecomposition.actualLU"].stop();
  timers["initializeSchurComplementSolver.LUDecomposition"].stop();
}
{'fraction_non_alphanumeric': 0.14936562860438293, 'fraction_numerical': 0.007991431866864393, 'mean_word_length': 3.280324400564175, 'pattern_counts': {'":': 0, '<': 128, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '145302', 'n_tokens_mistral': 4454, 'n_tokens_neox': 4261, 'n_words': 821}