oomph::CRDoubleMatrixHelpers Namespace Reference

Namespace for helper functions for CRDoubleMatrices. More...

Functions
void	create_uniformly_distributed_matrix (const unsigned &nrow, const unsigned &ncol, const OomphCommunicator *const comm_pt, const Vector< double > &values, const Vector< int > &column_indices, const Vector< int > &row_start, CRDoubleMatrix &matrix_out)
double	inf_norm (const DenseMatrix< CRDoubleMatrix * > &matrix_pt)
	Compute infinity (maximum) norm of sub blocks as if it was one matrix. More...
double	gershgorin_eigenvalue_estimate (const DenseMatrix< CRDoubleMatrix * > &matrix_pt)
void	concatenate (const DenseMatrix< CRDoubleMatrix * > &matrix_pt, CRDoubleMatrix &result_matrix)
void	concatenate_without_communication (const Vector< LinearAlgebraDistribution * > &row_distribution_pt, const Vector< LinearAlgebraDistribution * > &col_distribution_pt, const DenseMatrix< CRDoubleMatrix * > &matrix_pt, CRDoubleMatrix &result_matrix)
void	concatenate_without_communication (const Vector< LinearAlgebraDistribution * > &block_distribution_pt, const DenseMatrix< CRDoubleMatrix * > &matrix_pt, CRDoubleMatrix &result_matrix)
void	deep_copy (const CRDoubleMatrix *const in_matrix_pt, CRDoubleMatrix &out_matrix)
	Create a deep copy of the matrix pointed to by in_matrix_pt. More...

Detailed Description

Namespace for helper functions for CRDoubleMatrices.

Function Documentation

concatenate()

void oomph::CRDoubleMatrixHelpers::concatenate	(	const DenseMatrix< CRDoubleMatrix * > &	matrix_pt,
		CRDoubleMatrix &	result_matrix
	)

Concatenate CRDoubleMatrix matrices. The in matrices are concatenated such that the block structure of the in matrices are preserved in the result matrix. Communication between processors is required. If the block structure of the sub matrices does not need to be preserved, consider using CRDoubleMatrixHelpers::concatenate_without_communication(...).

The matrix manipulation functions CRDoubleMatrixHelpers::concatenate(...) and CRDoubleMatrixHelpers::concatenate_without_communication(...) are analogous to the Vector manipulation functions DoubleVectorHelpers::concatenate(...) and DoubleVectorHelpers::concatenate_without_communication(...). Please look at the DoubleVector functions for an illustration of the differences between concatenate(...) and concatenate_without_communication(...).

Distribution of the result matrix: If the result matrix does not have a distribution built, then it will be given a uniform row distribution. Otherwise we use the existing distribution. This gives the user the ability to define their own distribution, or save computing power if a distribution has been pre-built.

NOTE: ALL the matrices pointed to by matrix_pt has to be built. This is not the case with concatenate_without_communication(...)

    {
      // The number of block rows and block columns.
      unsigned matrix_nrow = matrix_pt.nrow();
      unsigned matrix_ncol = matrix_pt.ncol();
 
      // PARANOID checks involving only the in matrices.
#ifdef PARANOID
      // Are there matrices to concatenate?
      if (matrix_nrow == 0)
      {
        std::ostringstream error_message;
        error_message << "There are no matrices to concatenate.\n";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
 
      // Does this matrix need concatenating?
      if ((matrix_nrow == 1) && (matrix_ncol == 1))
      {
        std::ostringstream warning_message;
        warning_message << "There is only one matrix to concatenate...\n"
                        << "This does not require concatenating...\n";
        OomphLibWarning(warning_message.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
      }
 
      // Are all sub matrices built?
      for (unsigned block_row_i = 0; block_row_i < matrix_nrow; block_row_i++)
      {
        for (unsigned block_col_i = 0; block_col_i < matrix_ncol; block_col_i++)
        {
          if (!(matrix_pt(block_row_i, block_col_i)->built()))
          {
            std::ostringstream error_message;
            error_message << "The sub matrix (" << block_row_i << ","
                          << block_col_i << ")\n"
                          << "is not built. \n";
            throw OomphLibError(error_message.str(),
                                OOMPH_CURRENT_FUNCTION,
                                OOMPH_EXCEPTION_LOCATION);
          }
        }
      }
 
      // Do all dimensions of sub matrices "make sense"?
      // Compare the number of rows of each block matrix in a block row.
      for (unsigned block_row_i = 0; block_row_i < matrix_nrow; block_row_i++)
      {
        // Use the first column to compare against the rest.
        unsigned long current_block_nrow = matrix_pt(block_row_i, 0)->nrow();
 
        // Compare against columns 1 to matrix_ncol - 1
        for (unsigned block_col_i = 1; block_col_i < matrix_ncol; block_col_i++)
        {
          // Get the nrow for this sub block.
          unsigned long subblock_nrow =
            matrix_pt(block_row_i, block_col_i)->nrow();
 
          if (current_block_nrow != subblock_nrow)
          {
            std::ostringstream error_message;
            error_message << "The sub matrix (" << block_row_i << ","
                          << block_col_i << ")\n"
                          << "requires nrow = " << current_block_nrow
                          << ", but has nrow = " << subblock_nrow << ".\n";
            throw OomphLibError(error_message.str(),
                                OOMPH_CURRENT_FUNCTION,
                                OOMPH_EXCEPTION_LOCATION);
          }
        }
      }
 
      // Compare the number of columns of each block matrix in a block column.
      for (unsigned block_col_i = 0; block_col_i < matrix_ncol; block_col_i++)
      {
        // Use the first row to compare against the rest.
        unsigned long current_block_ncol = matrix_pt(0, block_col_i)->ncol();
 
        // Compare against rows 1 to matrix_nrow - 1
        for (unsigned block_row_i = 1; block_row_i < matrix_nrow; block_row_i++)
        {
          // Get the ncol for this sub block.
          unsigned long subblock_ncol =
            matrix_pt(block_row_i, block_col_i)->ncol();
 
          if (current_block_ncol != subblock_ncol)
          {
            std::ostringstream error_message;
            error_message << "The sub matrix (" << block_row_i << ","
                          << block_col_i << ")\n"
                          << "requires ncol = " << current_block_ncol
                          << ", but has ncol = " << subblock_ncol << ".\n";
            throw OomphLibError(error_message.str(),
                                OOMPH_CURRENT_FUNCTION,
                                OOMPH_EXCEPTION_LOCATION);
          }
        }
      }
#endif
 
      // The communicator pointer from block (0,0)
      const OomphCommunicator* const comm_pt =
        matrix_pt(0, 0)->distribution_pt()->communicator_pt();
 
      // Check if the block (0,0) is distributed or not.
      bool distributed = matrix_pt(0, 0)->distributed();
 
      // If the result matrix does not have a distribution, we create a uniform
      // distribution.
      if (!result_matrix.distribution_pt()->built())
      {
        // Sum of sub matrix nrow. We use the first column.
        unsigned tmp_nrow = 0;
        for (unsigned block_row_i = 0; block_row_i < matrix_nrow; block_row_i++)
        {
          tmp_nrow += matrix_pt(block_row_i, 0)->nrow();
        }
 
        LinearAlgebraDistribution tmp_distribution(
          comm_pt, tmp_nrow, distributed);
 
        result_matrix.build(&tmp_distribution);
      }
      else
      // A distribution is supplied for the result matrix.
      {
#ifdef PARANOID
        // Check if the sum of the nrow from the sub matrices is the same as the
        // the nrow from the result matrix.
 
        // Sum of sub matrix nrow. We use the first column.
        unsigned tmp_nrow = 0;
        for (unsigned block_row_i = 0; block_row_i < matrix_nrow; block_row_i++)
        {
          tmp_nrow += matrix_pt(block_row_i, 0)->nrow();
        }
 
        if (tmp_nrow != result_matrix.nrow())
        {
          std::ostringstream error_message;
          error_message << "The total number of rows from the matrices to\n"
                        << "concatenate does not match the nrow from the\n"
                        << "result matrix\n";
          throw OomphLibError(error_message.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
        }
#endif
      }
 
#ifdef PARANOID
 
      // Are all the communicators the same?
      // Compare the communicator for sub matrices (against the result matrix).
      {
        const OomphCommunicator communicator =
          *(result_matrix.distribution_pt()->communicator_pt());
 
        // Are all communicator pointers the same?
        for (unsigned block_row_i = 0; block_row_i < matrix_nrow; block_row_i++)
        {
          for (unsigned block_col_i = 0; block_col_i < matrix_ncol;
               block_col_i++)
          {
            const OomphCommunicator another_communicator =
              *(matrix_pt(block_row_i, block_col_i)
                  ->distribution_pt()
                  ->communicator_pt());
 
            if (!(communicator == another_communicator))
            {
              std::ostringstream error_message;
              error_message << "The OomphCommunicator of the sub matrix ("
                            << block_row_i << "," << block_col_i << ")\n"
                            << "does not have the same communicator as the "
                               "result matrix. \n";
              throw OomphLibError(error_message.str(),
                                  OOMPH_CURRENT_FUNCTION,
                                  OOMPH_EXCEPTION_LOCATION);
            }
          }
        }
      }
 
      // Are all the distributed boolean the same? This only applies if we have
      // more than one processor. If there is only one processor, then it does
      // not matter if it is distributed or not - they are conceptually the
      // same.
      if (comm_pt->nproc() != 1)
      {
        // Compare distributed for sub matrices (against the result matrix).
        const bool res_distributed = result_matrix.distributed();
 
        // Loop over all sub blocks.
        for (unsigned block_row_i = 0; block_row_i < matrix_nrow; block_row_i++)
        {
          for (unsigned block_col_i = 0; block_col_i < matrix_ncol;
               block_col_i++)
          {
            const bool another_distributed =
              matrix_pt(block_row_i, block_col_i)->distributed();
 
            if (res_distributed != another_distributed)
            {
              std::ostringstream error_message;
              error_message << "The distributed boolean of the sub matrix ("
                            << block_row_i << "," << block_col_i << ")\n"
                            << "is not the same as the result matrix. \n";
              throw OomphLibError(error_message.str(),
                                  OOMPH_CURRENT_FUNCTION,
                                  OOMPH_EXCEPTION_LOCATION);
            }
          }
        }
      }
#endif
 
 
      // Get the number of columns up to each block for offset
      // in calculating the result column indices.
      // Since the number of columns in each block column is the same,
      // we only loop through the first block row (row zero).
      Vector<unsigned long> sum_of_ncol_up_to_block(matrix_ncol);
 
      // Also compute the total number of columns to build the resulting matrix.
      unsigned long res_ncol = 0;
 
      for (unsigned block_col_i = 0; block_col_i < matrix_ncol; block_col_i++)
      {
        sum_of_ncol_up_to_block[block_col_i] = res_ncol;
        res_ncol += matrix_pt(0, block_col_i)->ncol();
      }
 
      // We begin the process of extracting and ordering the values,
      // column_indices and row_start of all the sub blocks.
      if ((comm_pt->nproc() == 1) || !distributed)
      // Serial version of the code.
      {
        // Get the total number of non zero entries so we can reserve storage
        // for the values and column_indices vectors.
        unsigned long res_nnz = 0;
        for (unsigned block_row_i = 0; block_row_i < matrix_nrow; block_row_i++)
        {
          for (unsigned block_col_i = 0; block_col_i < matrix_ncol;
               block_col_i++)
          {
            res_nnz += matrix_pt(block_row_i, block_col_i)->nnz();
          }
        }
 
        // Declare the vectors required to build a CRDoubleMatrix
        Vector<double> res_values;
        Vector<int> res_column_indices;
        Vector<int> res_row_start;
 
        // Reserve space for the vectors.
        res_values.reserve(res_nnz);
        res_column_indices.reserve(res_nnz);
        res_row_start.reserve(result_matrix.nrow() + 1);
 
        // Now we fill in the data.
 
        // Running sum of nnz per row.
        int nnz_running_sum = 0;
 
        // Loop through the block rows.
        for (unsigned block_row_i = 0; block_row_i < matrix_nrow; block_row_i++)
        {
          // Get the number of rows in this block row, from the first block.
          unsigned long block_row_nrow = matrix_pt(block_row_i, 0)->nrow();
 
          // Loop through the number of rows in this block row
          for (unsigned row_i = 0; row_i < block_row_nrow; row_i++)
          {
            // The row start is the nnz at the start of the row.
            res_row_start.push_back(nnz_running_sum);
 
            // Loop through the block columns
            for (unsigned block_col_i = 0; block_col_i < matrix_ncol;
                 block_col_i++)
            {
              // Get the current block.
              CRDoubleMatrix* current_block_pt =
                matrix_pt(block_row_i, block_col_i);
 
              // Get the values, column_indices and row_start for this block.
              double* current_block_values = current_block_pt->value();
              int* current_block_column_indices =
                current_block_pt->column_index();
              int* current_block_row_start = current_block_pt->row_start();
 
              for (int val_i = current_block_row_start[row_i];
                   val_i < current_block_row_start[row_i + 1];
                   val_i++)
              {
                res_values.push_back(current_block_values[val_i]);
                res_column_indices.push_back(
                  current_block_column_indices[val_i] +
                  sum_of_ncol_up_to_block[block_col_i]);
              }
 
              // Update the running sum of nnz per row
              nnz_running_sum += current_block_row_start[row_i + 1] -
                                 current_block_row_start[row_i];
            } // for block cols
          } // for rows
        } // for block rows
 
        // Fill in the last row start
        res_row_start.push_back(res_nnz);
 
        // Build the matrix
        result_matrix.build(
          res_ncol, res_values, res_column_indices, res_row_start);
      }
      // Otherwise we are dealing with a distributed matrix.
      else
      {
#ifdef OOMPH_HAS_MPI
 
        // Flag to enable timing. This is for debugging
        // and/or testing purposes only.
        bool enable_timing = false;
 
        // Get the number of processors
        unsigned nproc = comm_pt->nproc();
 
        // My rank
        unsigned my_rank = comm_pt->my_rank();
 
        // Storage for the data (per processor) to send.
        Vector<Vector<unsigned>> column_indices_to_send(nproc);
        Vector<Vector<double>> values_to_send(nproc);
 
        // The sum of the nrow for the sub blocks (so far). This is used as an
        // offset to calculate the global equation number in the result matrix.
        unsigned long sum_of_block_nrow = 0;
 
        double t_prep_data_start;
        if (enable_timing)
        {
          t_prep_data_start = TimingHelpers::timer();
        }
 
        // Get the pointer to the result distribution, for convenience...
        LinearAlgebraDistribution* res_distribution_pt =
          result_matrix.distribution_pt();
 
        // loop over the sub blocks to calculate the global_eqn, get the values
        // and column indices.
        for (unsigned block_row_i = 0; block_row_i < matrix_nrow; block_row_i++)
        {
          // Get the number of local rows in this block_row from the first
          // block.
          unsigned current_block_nrow_local =
            matrix_pt(block_row_i, 0)->nrow_local();
 
          // Get the first_row for this block_row
          unsigned current_block_row_first_row =
            matrix_pt(block_row_i, 0)->first_row();
 
          // Loop through the number of local rows
          for (unsigned sub_local_eqn = 0;
               sub_local_eqn < current_block_nrow_local;
               sub_local_eqn++)
          {
            // Calculate the corresponding (res_global_eqn) equation number
            // for this local row number in this block.
            unsigned long res_global_eqn =
              sub_local_eqn + current_block_row_first_row + sum_of_block_nrow;
 
            // Get the processor that this global row belongs to.
            // The rank_of_global_row(...) function loops through all the
            // processors and does two unsigned comparisons. Since we have to do
            // this for every row, it may be better to store a list mapping for
            // very large number of processors.
            unsigned res_p =
              res_distribution_pt->rank_of_global_row(res_global_eqn);
 
            // With the res_p, we get the res_first_row to
            // work out the res_local_eqn
            unsigned res_first_row = res_distribution_pt->first_row(res_p);
            unsigned res_local_eqn = res_global_eqn - res_first_row;
 
            // Loop through the block columns, calculate the nnz. This is used
            // to reserve space for the value and column_indices Vectors.
            unsigned long current_row_nnz = 0;
            for (unsigned block_col_i = 0; block_col_i < matrix_ncol;
                 block_col_i++)
            {
              // Get the row_start
              int* current_block_row_start =
                matrix_pt(block_row_i, block_col_i)->row_start();
 
              // Update the nnz for this row.
              current_row_nnz += current_block_row_start[sub_local_eqn + 1] -
                                 current_block_row_start[sub_local_eqn];
            } // for block column, get nnz.
 
            // Reserve space for efficiency.
            // unsigned capacity_in_res_p_vec
            //  = column_indices_to_send[res_p].capacity();
 
            // Reserve memory for nnz+2, since we need to store the
            // res_local_eqn and nnz as well as the data (values/column
            // indices). Note: The two reserve functions are called per row. If
            // the matrix is very sparse (just a few elements per row), it will
            // be more efficient to not reserve and let the STL vector handle
            // this. On average, this implementation is more efficient.
            // column_indices_to_send[res_p].reserve(capacity_in_res_p_vec
            //    + current_row_nnz+2);
            // values_to_send[res_p].reserve(capacity_in_res_p_vec
            //    + current_row_nnz+2);
 
            // Push back the res_local_eqn and nnz
            column_indices_to_send[res_p].push_back(res_local_eqn);
            column_indices_to_send[res_p].push_back(current_row_nnz);
            values_to_send[res_p].push_back(res_local_eqn);
            values_to_send[res_p].push_back(current_row_nnz);
 
            // Loop through the block columns again and get the values
            // and column_indices
            for (unsigned block_col_i = 0; block_col_i < matrix_ncol;
                 block_col_i++)
            {
              // Cache the pointer to the current block for convenience.
              CRDoubleMatrix* current_block_pt =
                matrix_pt(block_row_i, block_col_i);
 
              // Values, column indices and row_start for the current block.
              double* current_block_values = current_block_pt->value();
              int* current_block_column_indices =
                current_block_pt->column_index();
              int* current_block_row_start = current_block_pt->row_start();
 
              // Loop though the values and column_indices
              for (int val_i = current_block_row_start[sub_local_eqn];
                   val_i < current_block_row_start[sub_local_eqn + 1];
                   val_i++)
              {
                // Push back the value.
                values_to_send[res_p].push_back(current_block_values[val_i]);
 
                // Push back the (offset) column index.
                column_indices_to_send[res_p].push_back(
                  current_block_column_indices[val_i] +
                  sum_of_ncol_up_to_block[block_col_i]);
              } // for block columns
            } // for block column, get values and column_indices.
          } // for sub_local_eqn
 
          // update the sum_of_block_nrow
          sum_of_block_nrow += matrix_pt(block_row_i, 0)->nrow();
 
        } // for block_row
 
        if (enable_timing)
        {
          double t_prep_data_finish = TimingHelpers::timer();
          double t_prep_data_time = t_prep_data_finish - t_prep_data_start;
          oomph_info << "Time for prep data: " << t_prep_data_time << std::endl;
        }
 
        // Prepare to send data!
 
        // Storage for the number of data to be sent to each processor.
        Vector<int> send_n(nproc, 0);
 
        // Storage for all the values/column indices to be sent
        // to each processor.
        Vector<double> send_values_data;
        Vector<unsigned> send_column_indices_data;
 
        // Storage location within send_values_data
        // (and send_column_indices_data) for data to be sent to each processor.
        Vector<int> send_displacement(nproc, 0);
 
        double t_total_ndata_start;
        if (enable_timing) t_total_ndata_start = TimingHelpers::timer();
 
        // Get the total amount of data which needs to be sent, so we can
        // reserve space for it.
        unsigned total_ndata = 0;
        for (unsigned rank = 0; rank < nproc; rank++)
        {
          if (rank != my_rank)
          {
            total_ndata += values_to_send[rank].size();
          }
        }
 
        if (enable_timing)
        {
          double t_total_ndata_finish = TimingHelpers::timer();
          double t_total_ndata_time =
            t_total_ndata_finish - t_total_ndata_start;
          oomph_info << "Time for total_ndata: " << t_total_ndata_time
                     << std::endl;
        }
 
        double t_flat_pack_start;
        if (enable_timing) t_flat_pack_start = TimingHelpers::timer();
 
        // Now we don't have to re-allocate data/memory when push_back is
        // called. Nb. Using push_back without reserving memory may cause
        // multiple re-allocation behind the scenes, this is expensive.
        send_values_data.reserve(total_ndata);
        send_column_indices_data.reserve(total_ndata);
 
        // Loop over all the processors to "flat pack" the data for sending.
        for (unsigned rank = 0; rank < nproc; rank++)
        {
          // Set the offset for the current processor
          // This only has to be done once for both values and column indices.
          send_displacement[rank] = send_values_data.size();
 
          // Don't bother to do anything if
          // the processor in the loop is the current processor.
          if (rank != my_rank)
          {
            // Put the values into the send data vector.
            // n_data is the same for both values and column indices.
            unsigned n_data = values_to_send[rank].size();
            for (unsigned j = 0; j < n_data; j++)
            {
              send_values_data.push_back(values_to_send[rank][j]);
              send_column_indices_data.push_back(
                column_indices_to_send[rank][j]);
            } // for
          } // if rank != my_rank
 
          // Find the number of data to be added to the vector.
          // send_n is the same for both values and column indices.
          send_n[rank] = send_values_data.size() - send_displacement[rank];
        } // loop over processors
 
        if (enable_timing)
        {
          double t_flat_pack_finish = TimingHelpers::timer();
          double t_flat_pack_time = t_flat_pack_finish - t_flat_pack_start;
          oomph_info << "t_flat_pack_time: " << t_flat_pack_time << std::endl;
        }
 
        double t_sendn_start;
        if (enable_timing) t_sendn_start = TimingHelpers::timer();
 
        // Strorage for the number of data to be received from each processor
        Vector<int> receive_n(nproc, 0);
 
        MPI_Alltoall(&send_n[0],
                     1,
                     MPI_INT,
                     &receive_n[0],
                     1,
                     MPI_INT,
                     comm_pt->mpi_comm());
 
        if (enable_timing)
        {
          double t_sendn_finish = TimingHelpers::timer();
          double t_sendn_time = t_sendn_finish - t_sendn_start;
          oomph_info << "t_sendn_time: " << t_sendn_time << std::endl;
        }
 
 
        // Prepare the data to be received
        // by working out the displacement from the received data
        // receive_displacement is the same for both values and column indices.
        Vector<int> receive_displacement(nproc, 0);
        int receive_data_count = 0;
        for (unsigned rank = 0; rank < nproc; rank++)
        {
          receive_displacement[rank] = receive_data_count;
          receive_data_count += receive_n[rank];
        }
 
        // Now resize the receive buffer for all data from all processors.
        // Make sure that it has a size of at least one.
        if (receive_data_count == 0)
        {
          receive_data_count++;
        }
        Vector<double> receive_values_data(receive_data_count);
        Vector<unsigned> receive_column_indices_data(receive_data_count);
 
        // Make sure that the send buffer has size at least one
        // so that we don't get a segmentation fault.
        if (send_values_data.size() == 0)
        {
          send_values_data.resize(1);
        }
 
        double t_send_data_start;
        if (enable_timing) t_send_data_start = TimingHelpers::timer();
 
        // Now send the data between all processors
        MPI_Alltoallv(&send_values_data[0],
                      &send_n[0],
                      &send_displacement[0],
                      MPI_DOUBLE,
                      &receive_values_data[0],
                      &receive_n[0],
                      &receive_displacement[0],
                      MPI_DOUBLE,
                      comm_pt->mpi_comm());
 
        // Now send the data between all processors
        MPI_Alltoallv(&send_column_indices_data[0],
                      &send_n[0],
                      &send_displacement[0],
                      MPI_UNSIGNED,
                      &receive_column_indices_data[0],
                      &receive_n[0],
                      &receive_displacement[0],
                      MPI_UNSIGNED,
                      comm_pt->mpi_comm());
 
        if (enable_timing)
        {
          double t_send_data_finish = TimingHelpers::timer();
          double t_send_data_time = t_send_data_finish - t_send_data_start;
          oomph_info << "t_send_data_time: " << t_send_data_time << std::endl;
        }
 
        // All the rows for this processor are stored in:
        // from other processors:
        // receive_column_indices_data and receive_values_data
        // from this processor:
        // column_indices_to_send[my_rank] and values_to_send[my_rank]
        //
        // They are in some order determined by the distribution.
        // We need to re-arrange them. To do this, we do some pre-processing.
 
        // nrow_local for this processor.
        unsigned res_nrow_local = res_distribution_pt->nrow_local();
 
        // Per row, store:
        // 1) where this row came from, 0 - this proc, 1 - other procs.
        // 2) the nnz,
        // 3) the offset - where the values/columns in the receive data vectors
        //                 begins. This is different from the offset of where
        //                 the data from a certain processor starts.
        Vector<Vector<unsigned>> value_column_locations(res_nrow_local,
                                                        Vector<unsigned>(3, 0));
 
        // Store the local nnz so we can reserve space for
        // the values and column indices.
        unsigned long res_nnz_local = 0;
 
        double t_locations_start;
        if (enable_timing) t_locations_start = TimingHelpers::timer();
 
        // Loop through the data currently on this processor.
        unsigned location_i = 0;
        unsigned my_column_indices_to_send_size =
          column_indices_to_send[my_rank].size();
        while (location_i < my_column_indices_to_send_size)
        {
          unsigned current_local_eqn =
            column_indices_to_send[my_rank][location_i++];
 
          unsigned current_nnz = column_indices_to_send[my_rank][location_i++];
 
          // No need to fill [*][0] with 0 since it is already initialised to 0.
 
          // Store the nnz.
          value_column_locations[current_local_eqn][1] = current_nnz;
 
          // Also increment the res_local_nnz
          res_nnz_local += current_nnz;
 
          // Store the offset.
          value_column_locations[current_local_eqn][2] = location_i;
 
          // Update the location_i so it starts at the next row.
          location_i += current_nnz;
        }
 
        // Loop through the data from different processors.
 
        // Check to see if data has been received.
        bool data_has_been_received = false;
        unsigned send_rank = 0;
        while (send_rank < nproc)
        {
          if (receive_n[send_rank] > 0)
          {
            data_has_been_received = true;
            break;
          }
          send_rank++;
        }
 
        location_i = 0; // start at 0.
        if (data_has_been_received)
        {
          unsigned receive_column_indices_data_size =
            receive_column_indices_data.size();
          while (location_i < receive_column_indices_data_size)
          {
            unsigned current_local_eqn =
              receive_column_indices_data[location_i++];
            unsigned current_nnz = receive_column_indices_data[location_i++];
 
            // These comes from other processors.
            value_column_locations[current_local_eqn][0] = 1;
 
            // Store the nnz.
            value_column_locations[current_local_eqn][1] = current_nnz;
 
            // Also increment the res_local_nnz
            res_nnz_local += current_nnz;
 
            // Store the offset.
            value_column_locations[current_local_eqn][2] = location_i;
 
            // Update the location_i so it starts at the next row.
            location_i += current_nnz;
          }
        }
 
        if (enable_timing)
        {
          double t_locations_finish = TimingHelpers::timer();
          double t_locations_time = t_locations_finish - t_locations_start;
          oomph_info << "t_locations_time: " << t_locations_time << std::endl;
        }
 
        double t_fillvecs_start;
        if (enable_timing) t_fillvecs_start = TimingHelpers::timer();
 
        // Now loop through the locations and store the values
        // the column indices in the correct order.
        Vector<int> res_column_indices;
        Vector<double> res_values;
        Vector<int> res_row_start;
 
        res_column_indices.reserve(res_nnz_local);
        res_values.reserve(res_nnz_local);
        res_row_start.reserve(res_nrow_local + 1);
 
        // Running sum of nnz for the row_start. Must be int because
        // res_row_start is templated with int.
        int nnz_running_sum = 0;
 
        // Now insert the rows.
        for (unsigned local_row_i = 0; local_row_i < res_nrow_local;
             local_row_i++)
        {
          // Fill the res_row_start with the nnz so far.
          res_row_start.push_back(nnz_running_sum);
 
          bool data_is_from_other_proc =
            bool(value_column_locations[local_row_i][0]);
 
          unsigned row_i_nnz = value_column_locations[local_row_i][1];
 
          unsigned row_i_offset = value_column_locations[local_row_i][2];
 
          if (data_is_from_other_proc)
          {
            // Insert range [offset, offset+nnz) from
            // receive_column_indices_data and receive_values_data into
            // res_column_indices and res_values respectively.
            res_column_indices.insert(
              res_column_indices.end(),
              receive_column_indices_data.begin() + row_i_offset,
              receive_column_indices_data.begin() + row_i_offset + row_i_nnz);
 
            res_values.insert(res_values.end(),
                              receive_values_data.begin() + row_i_offset,
                              receive_values_data.begin() + row_i_offset +
                                row_i_nnz);
          }
          else
          {
            res_column_indices.insert(res_column_indices.end(),
                                      column_indices_to_send[my_rank].begin() +
                                        row_i_offset,
                                      column_indices_to_send[my_rank].begin() +
                                        row_i_offset + row_i_nnz);
 
            res_values.insert(res_values.end(),
                              values_to_send[my_rank].begin() + row_i_offset,
                              values_to_send[my_rank].begin() + row_i_offset +
                                row_i_nnz);
          }
 
          // Update the running sum of nnz
          nnz_running_sum += row_i_nnz;
        }
 
        // Insert the last row_start value
        res_row_start.push_back(res_nnz_local);
 
        if (enable_timing)
        {
          double t_fillvecs_finish = TimingHelpers::timer();
          double t_fillvecs_time = t_fillvecs_finish - t_fillvecs_start;
          oomph_info << "t_fillvecs_time: " << t_fillvecs_time << std::endl;
        }
 
        double t_buildres_start;
        if (enable_timing) t_buildres_start = TimingHelpers::timer();
 
        // build the matrix.
        result_matrix.build(
          res_ncol, res_values, res_column_indices, res_row_start);
 
        if (enable_timing)
        {
          double t_buildres_finish = TimingHelpers::timer();
          double t_buildres_time = t_buildres_finish - t_buildres_start;
          oomph_info << "t_buildres_time: " << t_buildres_time << std::endl;
        }
        //  */
#endif
      }
    }

References oomph::CRDoubleMatrix::build(), oomph::LinearAlgebraDistribution::built(), oomph::CRDoubleMatrix::column_index(), oomph::LinearAlgebraDistribution::communicator_pt(), oomph::DistributableLinearAlgebraObject::distributed(), oomph::DistributableLinearAlgebraObject::distribution_pt(), oomph::LinearAlgebraDistribution::first_row(), j, oomph::OomphCommunicator::my_rank(), oomph::DenseMatrix< T >::ncol(), oomph::OomphCommunicator::nproc(), oomph::DenseMatrix< T >::nrow(), oomph::CRDoubleMatrix::nrow(), oomph::LinearAlgebraDistribution::nrow_local(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, oomph::LinearAlgebraDistribution::rank_of_global_row(), oomph::CRDoubleMatrix::row_start(), oomph::TimingHelpers::timer(), and oomph::CRDoubleMatrix::value().

Referenced by main(), and oomph::SpaceTimeNavierStokesSubsidiaryPreconditioner::setup().

concatenate_without_communication() [1/2]

void oomph::CRDoubleMatrixHelpers::concatenate_without_communication	(	const Vector< LinearAlgebraDistribution * > &	block_distribution_pt,
		const DenseMatrix< CRDoubleMatrix * > &	matrix_pt,
		CRDoubleMatrix &	result_matrix
	)

Concatenate CRDoubleMatrix matrices. This calls the other concatenate_without_communication(...) function, passing block_distribution_pt as both the row_distribution_pt and col_distribution_pt. This should only be called for block square matrices.

    {
#ifdef PARANOID
      // The number of block rows and block columns.
      unsigned matrix_nrow = matrix_pt.nrow();
      unsigned matrix_ncol = matrix_pt.ncol();
 
      // Are there matrices to concatenate?
      if (matrix_nrow == 0)
      {
        std::ostringstream error_message;
        error_message << "There are no matrices to concatenate.\n";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
 
      // Ensure that the sub matrices is a square block matrix.
      if (matrix_nrow != matrix_ncol)
      {
        std::ostringstream error_message;
        error_message
          << "The number of block rows and block columns\n"
          << "must be the same. Otherwise, call the other\n"
          << "concatenate_without_communication function, passing in\n"
          << "a Vector of distributions describing how to permute the\n"
          << "columns.";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      concatenate_without_communication(
        block_distribution_pt, block_distribution_pt, matrix_pt, result_matrix);
    }

References concatenate_without_communication(), oomph::DenseMatrix< T >::ncol(), oomph::DenseMatrix< T >::nrow(), OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

concatenate_without_communication() [2/2]

void oomph::CRDoubleMatrixHelpers::concatenate_without_communication	(	const Vector< LinearAlgebraDistribution * > &	row_distribution_pt,
		const Vector< LinearAlgebraDistribution * > &	col_distribution_pt,
		const DenseMatrix< CRDoubleMatrix * > &	matrix_pt,
		CRDoubleMatrix &	result_matrix
	)

Concatenate CRDoubleMatrix matrices.

The Vector row_distribution_pt contains the LinearAlgebraDistribution of each block row. The Vector col_distribution_pt contains the LinearAlgebraDistribution of each block column. The DenseMatrix matrix_pt contains pointers to the CRDoubleMatrices to concatenate. The CRDoubleMatrix result_matrix is the result matrix.

The result matrix is a permutation of the sub matrices such that the data stays on the same processor when the result matrix is built, there is no communication between processors. Thus the block structure of the sub matrices are NOT preserved in the result matrix. The rows are block-permuted, defined by the concatenation of the distributions in row_distribution_pt. Similarly, the columns are block-permuted, defined by the concatenation of the distributions in col_distribution_pt. For more details on the block-permutation, see LinearAlgebraDistributionHelpers::concatenate(...).

If one wishes to preserve the block structure of the sub matrices in the result matrix, consider using CRDoubleMatrixHelpers::concatenate(...), which uses communication between processors to ensure that the block structure of the sub matrices are preserved.

The matrix manipulation functions CRDoubleMatrixHelpers::concatenate(...) and CRDoubleMatrixHelpers::concatenate_without_communication(...) are analogous to the Vector manipulation functions DoubleVectorHelpers::concatenate(...) and DoubleVectorHelpers::concatenate_without_communication(...). Please look at the DoubleVector functions for an illustration of the differences between concatenate(...) and concatenate_without_communication(...).

Distribution of the result matrix: If the result matrix does not have a distribution built, then it will be given a distribution built from the concatenation of the distributions from row_distribution_pt, see LinearAlgebraDistributionHelpers::concatenate(...) for more detail. Otherwise we use the existing distribution. If there is an existing distribution then it must be the same as the distribution from the concatenation of row distributions as described above. Why don't we always compute the distribution "on the fly"? Because a non-uniform distribution requires communication. All block preconditioner distributions are concatenations of the distributions of the individual blocks.

/////////////// END OF PARANOID TESTS ////////////////////////////////////

    {
      // The number of block rows and block columns.
      unsigned matrix_nrow = matrix_pt.nrow();
      unsigned matrix_ncol = matrix_pt.ncol();
 
      // PARANOID checks involving in matrices and block_distribution only.
      // PARANOID checks involving the result matrix will come later since
      // we have to create the result matrix distribution from the in
      // distribution if it does not already exist.
#ifdef PARANOID
 
      // Are there matrices to concatenate?
      if (matrix_nrow == 0 || matrix_ncol == 0)
      {
        std::ostringstream error_message;
        error_message << "There are no matrices to concatenate.\n";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
 
      // Does this matrix need concatenating?
      if ((matrix_nrow == 1) && (matrix_ncol == 1))
      {
        std::ostringstream warning_message;
        warning_message << "There is only one matrix to concatenate...\n"
                        << "This does not require concatenating...\n";
        OomphLibWarning(warning_message.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
      }
 
 
      // The distribution for each block row is stored in row_distribution_pt.
      // So the number of distributions in row_distribution_pt must be the
      // same as matrix_nrow.
      if (matrix_nrow != row_distribution_pt.size())
      {
        std::ostringstream error_message;
        error_message << "The number of row distributions must be the same as\n"
                      << "the number of block rows.";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
 
      // The number of distributions for the columns must match the number of
      // block columns.
      if (matrix_ncol != col_distribution_pt.size())
      {
        std::ostringstream error_message;
        error_message
          << "The number of column distributions must be the same as\n"
          << "the number of block columns.";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
 
      // Check that all pointers in row_distribution_pt is not null.
      for (unsigned block_row_i = 0; block_row_i < matrix_nrow; block_row_i++)
      {
        if (row_distribution_pt[block_row_i] == 0)
        {
          std::ostringstream error_message;
          error_message << "The row distribution pointer in position "
                        << block_row_i << " is null.\n";
          throw OomphLibError(error_message.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
        }
      }
 
      // Check that all pointers in row_distribution_pt is not null.
      for (unsigned block_col_i = 0; block_col_i < matrix_ncol; block_col_i++)
      {
        if (col_distribution_pt[block_col_i] == 0)
        {
          std::ostringstream error_message;
          error_message << "The column distribution pointer in position "
                        << block_col_i << " is null.\n";
          throw OomphLibError(error_message.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
        }
      }
 
      // Check that all distributions are built.
      // First the row distributions
      for (unsigned block_row_i = 0; block_row_i < matrix_nrow; block_row_i++)
      {
        if (!row_distribution_pt[block_row_i]->built())
        {
          std::ostringstream error_message;
          error_message << "The distribution pointer in position "
                        << block_row_i << " is not built.\n";
          throw OomphLibError(error_message.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
        }
      }
      // Now the column distributions
      for (unsigned block_col_i = 0; block_col_i < matrix_ncol; block_col_i++)
      {
        if (!col_distribution_pt[block_col_i]->built())
        {
          std::ostringstream error_message;
          error_message << "The distribution pointer in position "
                        << block_col_i << " is not built.\n";
          throw OomphLibError(error_message.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
        }
      }
 
      // Check that all communicators in row_distribution_pt are the same.
      const OomphCommunicator first_row_comm =
        *(row_distribution_pt[0]->communicator_pt());
 
      for (unsigned block_row_i = 1; block_row_i < matrix_nrow; block_row_i++)
      {
        const OomphCommunicator current_comm =
          *(row_distribution_pt[block_row_i]->communicator_pt());
 
        if (first_row_comm != current_comm)
        {
          std::ostringstream error_message;
          error_message
            << "The communicator from the row distribution in position "
            << block_row_i << " is not the same as the first "
            << "communicator from row_distribution_pt";
          throw OomphLibError(error_message.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
        }
      }
 
      // Check that all communicators in col_distribution_pt are the same as the
      // first row communicator from above.
      for (unsigned block_col_i = 0; block_col_i < matrix_ncol; block_col_i++)
      {
        const OomphCommunicator current_comm =
          *(col_distribution_pt[block_col_i]->communicator_pt());
 
        if (first_row_comm != current_comm)
        {
          std::ostringstream error_message;
          error_message
            << "The communicator from the col distribution in position "
            << block_col_i << " is not the same as the first "
            << "communicator from row_distribution_pt";
          throw OomphLibError(error_message.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
        }
      }
 
      // Are all sub matrices built? If the matrix_pt is not null, make sure
      // that it is built.
      for (unsigned block_row_i = 0; block_row_i < matrix_nrow; block_row_i++)
      {
        for (unsigned block_col_i = 0; block_col_i < matrix_ncol; block_col_i++)
        {
          if (matrix_pt(block_row_i, block_col_i) != 0 &&
              !(matrix_pt(block_row_i, block_col_i)->built()))
          {
            std::ostringstream error_message;
            error_message << "The sub matrix_pt(" << block_row_i << ","
                          << block_col_i << ")\n"
                          << "is not built.\n";
            throw OomphLibError(error_message.str(),
                                OOMPH_CURRENT_FUNCTION,
                                OOMPH_EXCEPTION_LOCATION);
          }
        }
      }
 
      // For the matrices which are built, do they have the same communicator as
      // the first communicator from row_distribution_pt?
      for (unsigned block_row_i = 0; block_row_i < matrix_nrow; block_row_i++)
      {
        for (unsigned block_col_i = 0; block_col_i < matrix_ncol; block_col_i++)
        {
          if (matrix_pt(block_row_i, block_col_i) != 0)
          {
            const OomphCommunicator current_comm =
              *(matrix_pt(block_row_i, block_col_i)
                  ->distribution_pt()
                  ->communicator_pt());
            if (first_row_comm != current_comm)
            {
              std::ostringstream error_message;
              error_message
                << "The sub matrix_pt(" << block_row_i << "," << block_col_i
                << ")\n"
                << "does not have the same communicator pointer as those in\n"
                << "(row|col)_distribution_pt.\n";
              throw OomphLibError(error_message.str(),
                                  OOMPH_CURRENT_FUNCTION,
                                  OOMPH_EXCEPTION_LOCATION);
            }
          }
        }
      }
 
      // Do all dimensions of sub matrices "make sense"?
      // Compare the number of rows of each block matrix in a block row.
      for (unsigned block_row_i = 0; block_row_i < matrix_nrow; block_row_i++)
      {
        // Use the first column to compare against the rest.
        unsigned long current_block_nrow =
          row_distribution_pt[block_row_i]->nrow();
 
        // Compare against columns 0 to matrix_ncol - 1
        for (unsigned block_col_i = 0; block_col_i < matrix_ncol; block_col_i++)
        {
          // Perform the check if the matrix_pt is not null.
          if (matrix_pt(block_row_i, block_col_i) != 0)
          {
            // Get the nrow for this sub block.
            unsigned long subblock_nrow =
              matrix_pt(block_row_i, block_col_i)->nrow();
 
            if (current_block_nrow != subblock_nrow)
            {
              std::ostringstream error_message;
              error_message << "The sub matrix (" << block_row_i << ","
                            << block_col_i << ")\n"
                            << "requires nrow = " << current_block_nrow
                            << ", but has nrow = " << subblock_nrow << ".\n"
                            << "Either the row_distribution_pt is incorrect or "
                            << "the sub matrices are incorrect.\n";
              throw OomphLibError(error_message.str(),
                                  OOMPH_CURRENT_FUNCTION,
                                  OOMPH_EXCEPTION_LOCATION);
            }
          }
        }
      }
 
      // Compare the number of columns of each block matrix in a block column.
      for (unsigned block_col_i = 0; block_col_i < matrix_ncol; block_col_i++)
      {
        // Get the current block ncol from the linear algebra distribution.
        // Note that we assume that the dimensions are symmetrical.
        unsigned current_block_ncol = col_distribution_pt[block_col_i]->nrow();
 
        for (unsigned block_row_i = 0; block_row_i < matrix_nrow; block_row_i++)
        {
          if (matrix_pt(block_row_i, block_col_i) != 0)
          {
            // Get the ncol for this sub block.
            unsigned subblock_ncol =
              matrix_pt(block_row_i, block_col_i)->ncol();
 
            if (current_block_ncol != subblock_ncol)
            {
              std::ostringstream error_message;
              error_message << "The sub matrix (" << block_row_i << ","
                            << block_col_i << ")\n"
                            << "requires ncol = " << current_block_ncol
                            << ", but has ncol = " << subblock_ncol << ".\n"
                            << "Either the col_distribution_pt is incorrect or "
                            << "the sub matrices are incorrect.\n";
              throw OomphLibError(error_message.str(),
                                  OOMPH_CURRENT_FUNCTION,
                                  OOMPH_EXCEPTION_LOCATION);
            }
          }
        }
      }
 
      // Ensure that the distributions for all sub matrices in the same block
      // row are the same. This is because we permute the row across several
      // matrices.
 
      // Loop through each block row.
      for (unsigned block_row_i = 0; block_row_i < matrix_nrow; block_row_i++)
      {
        // Get the distribution from the first block in this row.
        LinearAlgebraDistribution* block_row_distribution_pt =
          row_distribution_pt[block_row_i];
 
        // Loop through the block columns
        for (unsigned block_col_i = 0; block_col_i < matrix_ncol; block_col_i++)
        {
          if (matrix_pt(block_row_i, block_col_i) != 0)
          {
            // Get the distribution for this block.
            LinearAlgebraDistribution* current_block_distribution_pt =
              matrix_pt(block_row_i, block_col_i)->distribution_pt();
 
            // Ensure that the in matrices is a square block matrix.
            if ((*block_row_distribution_pt) !=
                (*current_block_distribution_pt))
            {
              std::ostringstream error_message;
              error_message
                << "Sub block(" << block_row_i << "," << block_col_i << ")"
                << "does not have the same distributoin as the first"
                << "block in this block row.\n"
                << "All distributions on a block row must be the same"
                << "for this function to concatenate matrices.\n";
              throw OomphLibError(error_message.str(),
                                  OOMPH_CURRENT_FUNCTION,
                                  OOMPH_EXCEPTION_LOCATION);
            }
          }
        }
      }
#endif
 
      // The communicator pointer from the first row_distribution_pt
      const OomphCommunicator* const comm_pt =
        row_distribution_pt[0]->communicator_pt();
 
      // Renamed for so it makes more sense.
      unsigned nblock_row = matrix_nrow;
 
      // If the result matrix does not have a distribution, then we concatenate
      // the sub matrix distributions.
      if (!result_matrix.distribution_pt()->built())
      {
        // The result distribution
        LinearAlgebraDistribution tmp_distribution;
        LinearAlgebraDistributionHelpers::concatenate(row_distribution_pt,
                                                      tmp_distribution);
 
        result_matrix.build(&tmp_distribution);
      }
      else
      // A distribution is supplied for the result matrix.
      {
#ifdef PARANOID
        // Check that the result distribution is a concatenation of the
        // distributions of the sub matrices.
 
        LinearAlgebraDistribution wanted_distribution;
 
        LinearAlgebraDistributionHelpers::concatenate(row_distribution_pt,
                                                      wanted_distribution);
 
        if (*(result_matrix.distribution_pt()) != wanted_distribution)
        {
          std::ostringstream error_message;
          error_message
            << "The result distribution is not correct.\n"
            << "Please call the function without a result\n"
            << "distribution (clear the result matrix) or check the\n"
            << "distribution of the result matrix.\n"
            << "The result distribution must be the same as the one \n"
            << "created by\n"
            << "LinearAlgebraDistributionHelpers::concatenate(...)";
          throw OomphLibError(error_message.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
        }
#endif
      }
 
      // The rest of the paranoid checks.
#ifdef PARANOID
 
      // Make sure that the communicator from the result matrix is the same as
      // all the others. This test is redundant if this function created the
      // result matrix distribution, since then it is guaranteed that the
      // communicators are the same.
      {
        // Communicator from the result matrix.
        const OomphCommunicator res_comm =
          *(result_matrix.distribution_pt()->communicator_pt());
 
        // Is the result communicator pointer the same as the others?
        // Since we have already tested the others, we only need to compare
        // against one of them. Say the first communicator from
        // row_distribution_pt.
        const OomphCommunicator first_comm =
          *(row_distribution_pt[0]->communicator_pt());
 
        if (res_comm != first_comm)
        {
          std::ostringstream error_message;
          error_message << "The OomphCommunicator of the result matrix is not "
                           "the same as the "
                        << "others!";
          throw OomphLibError(error_message.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
        }
      }
 
      // Are all the distributed boolean the same? This only applies if we have
      // more than one processor. If there is only one processor, then it does
      // not matter if it is distributed or not - they are conceptually the
      // same.
      if (comm_pt->nproc() != 1)
      {
        // Compare distributed for sub matrices (against the result matrix).
        const bool res_distributed = result_matrix.distributed();
 
        // Loop over all sub blocks.
        for (unsigned block_row_i = 0; block_row_i < matrix_nrow; block_row_i++)
        {
          for (unsigned block_col_i = 0; block_col_i < matrix_ncol;
               block_col_i++)
          {
            if (matrix_pt(block_row_i, block_col_i) != 0)
            {
              const bool another_distributed =
                matrix_pt(block_row_i, block_col_i)->distributed();
 
              if (res_distributed != another_distributed)
              {
                std::ostringstream error_message;
                error_message << "The distributed boolean of the sub matrix ("
                              << block_row_i << "," << block_col_i << ")\n"
                              << "is not the same as the result matrix. \n";
                throw OomphLibError(error_message.str(),
                                    OOMPH_CURRENT_FUNCTION,
                                    OOMPH_EXCEPTION_LOCATION);
              }
            }
          }
        }
 
        // Do this test for row_distribution_pt
        const bool first_row_distribution_distributed =
          row_distribution_pt[0]->distributed();
 
        for (unsigned block_row_i = 1; block_row_i < matrix_nrow; block_row_i++)
        {
          const bool another_distributed =
            row_distribution_pt[block_row_i]->distributed();
 
          if (first_row_distribution_distributed != another_distributed)
          {
            std::ostringstream error_message;
            error_message
              << "The distributed boolean of row_distribution_pt["
              << block_row_i << "]\n"
              << "is not the same as the one from row_distribution_pt[0]. \n";
            throw OomphLibError(error_message.str(),
                                OOMPH_CURRENT_FUNCTION,
                                OOMPH_EXCEPTION_LOCATION);
          }
        }
 
        // Repeat for col_distribution_pt
        for (unsigned block_col_i = 0; block_col_i < matrix_ncol; block_col_i++)
        {
          const bool another_distributed =
            col_distribution_pt[block_col_i]->distributed();
 
          if (first_row_distribution_distributed != another_distributed)
          {
            std::ostringstream error_message;
            error_message
              << "The distributed boolean of col_distribution_pt["
              << block_col_i << "]\n"
              << "is not the same as the one from row_distribution_pt[0]. \n";
            throw OomphLibError(error_message.str(),
                                OOMPH_CURRENT_FUNCTION,
                                OOMPH_EXCEPTION_LOCATION);
          }
        }
      }
#endif
 
 
      // The number of processors.
      unsigned nproc = comm_pt->nproc();
 
      // Cache the result distribution pointer for convenience.
      LinearAlgebraDistribution* res_distribution_pt =
        result_matrix.distribution_pt();
 
      // nrow_local for the result matrix
      unsigned res_nrow_local = res_distribution_pt->nrow_local();
 
      // renamed for readability.
      unsigned nblock_col = matrix_ncol;
 
      // construct the block offset
      //  DenseMatrix<unsigned> col_offset(nproc,nblock_col,0);
      std::vector<std::vector<unsigned>> col_offset(
        nproc, std::vector<unsigned>(nblock_col));
      unsigned off = 0;
      for (unsigned proc_i = 0; proc_i < nproc; proc_i++)
      {
        for (unsigned block_i = 0; block_i < nblock_col; block_i++)
        {
          col_offset[proc_i][block_i] = off;
          off += col_distribution_pt[block_i]->nrow_local(proc_i);
        }
      }
 
      // Do some pre-processing for the processor number a global row number is
      // on. This is required when permuting the column entries.
      // We need to do this for each distribution, so we have a vector of
      // vectors. First index corresponds to the distribution, the second is
      // the processor number.
      std::vector<std::vector<unsigned>> p_for_rows(nblock_col,
                                                    std::vector<unsigned>());
      // initialise 2D vector
      for (unsigned blocki = 0; blocki < nblock_col; blocki++)
      {
        int blockinrow = col_distribution_pt[blocki]->nrow();
        p_for_rows[blocki].resize(blockinrow);
        // FOR each global index in the block, work out the corresponding proc.
        for (int rowi = 0; rowi < blockinrow; rowi++)
        {
          unsigned p = 0;
          int b_first_row = col_distribution_pt[blocki]->first_row(p);
          int b_nrow_local = col_distribution_pt[blocki]->nrow_local(p);
 
          while (rowi < b_first_row || rowi >= b_nrow_local + b_first_row)
          {
            p++;
            b_first_row = col_distribution_pt[blocki]->first_row(p);
            b_nrow_local = col_distribution_pt[blocki]->nrow_local(p);
          }
          p_for_rows[blocki][rowi] = p;
        }
      }
 
      // determine nnz of all blocks on this processor only.
      // This is used to create storage space.
      unsigned long res_nnz = 0;
      for (unsigned row_i = 0; row_i < nblock_row; row_i++)
      {
        for (unsigned col_i = 0; col_i < nblock_col; col_i++)
        {
          if (matrix_pt(row_i, col_i) != 0)
          {
            res_nnz += matrix_pt(row_i, col_i)->nnz();
          }
        }
      }
 
      // My rank
      //    unsigned my_rank = comm_pt->my_rank();
      //    my_rank = my_rank;
 
      // Turn the above into a string.
      //      std::ostringstream myrankstream;
      //      myrankstream << "THISDOESNOTHINGnp" << my_rank << std::endl;
      //      std::string myrankstring = myrankstream.str();
 
 
      // CALLGRIND_ZERO_STATS;
      // CALLGRIND_START_INSTRUMENTATION;
 
      // storage for the result matrix.
      int* res_row_start = new int[res_nrow_local + 1];
      int* res_column_index = new int[res_nnz];
      double* res_value = new double[res_nnz];
 
      // initialise the zero-th entry
      res_row_start[0] = 0;
 
      // loop over the block rows
      unsigned long res_i = 0; // index for the result matrix.
      unsigned long res_row_i = 0; // index for the row
      for (unsigned i = 0; i < nblock_row; i++)
      {
        // loop over the rows of the current block local rows.
        unsigned block_nrow = row_distribution_pt[i]->nrow_local();
        for (unsigned k = 0; k < block_nrow; k++)
        {
          // initialise res_row_start
          res_row_start[res_row_i + 1] = res_row_start[res_row_i];
 
          // Loop over the block columns
          for (unsigned j = 0; j < nblock_col; j++)
          {
            // if block(i,j) pointer is not null then
            if (matrix_pt(i, j) != 0)
            {
              // get pointers for the elements in the current block
              int* b_row_start = matrix_pt(i, j)->row_start();
              int* b_column_index = matrix_pt(i, j)->column_index();
              double* b_value = matrix_pt(i, j)->value();
 
              // memcpy( &dst[dstIdx], &src[srcIdx], numElementsToCopy * sizeof(
              // Element ) );
              // no ele to copy
              int numEleToCopy = b_row_start[k + 1] - b_row_start[k];
              memcpy(res_value + res_i,
                     b_value + b_row_start[k],
                     numEleToCopy * sizeof(double));
              // Loop through the current local row.
              for (int l = b_row_start[k]; l < b_row_start[k + 1]; l++)
              {
                // if b_column_index[l] was a row index, what processor
                // would it be on
                //            unsigned p = col_distribution_pt[j]
                //              ->rank_of_global_row_map(b_column_index[l]);
                unsigned p = p_for_rows[j][b_column_index[l]];
 
                int b_first_row = col_distribution_pt[j]->first_row(p);
                //            int b_nrow_local =
                //            col_distribution_pt[j]->nrow_local(p);
 
                //            while (b_column_index[l] < b_first_row ||
                //                   b_column_index[l] >=
                //                   b_nrow_local+b_first_row)
                //             {
                //              p++;
                //              b_first_row =
                //              col_distribution_pt[j]->first_row(p);
                //              b_nrow_local =
                //              col_distribution_pt[j]->nrow_local(p);
                //             }
 
                // determine the local equation number in the block j/processor
                // p "column block"
                int eqn = b_column_index[l] - b_first_row;
 
                // add to the result matrix
                //            res_value[res_i] = b_value[l];
                res_column_index[res_i] = col_offset[p][j] + eqn;
                res_row_start[res_row_i + 1]++;
                res_i++;
              }
            }
          }
 
          // increment the row pt
          res_row_i++;
        }
      }
      // CALLGRIND_STOP_INSTRUMENTATION;
      // CALLGRIND_DUMP_STATS_AT(myrankstring.c_str());
 
 
      // Get the number of columns of the result matrix.
      unsigned res_ncol = 0;
      for (unsigned block_col_i = 0; block_col_i < matrix_ncol; block_col_i++)
      {
        res_ncol += col_distribution_pt[block_col_i]->nrow();
      }
 
      // Build the result matrix.
      result_matrix.build_without_copy(
        res_ncol, res_nnz, res_value, res_column_index, res_row_start);
    }

References oomph::CRDoubleMatrix::build(), oomph::CRDoubleMatrix::build_without_copy(), oomph::LinearAlgebraDistribution::built(), oomph::LinearAlgebraDistribution::communicator_pt(), oomph::LinearAlgebraDistributionHelpers::concatenate(), oomph::DistributableLinearAlgebraObject::distributed(), oomph::DistributableLinearAlgebraObject::distribution_pt(), i, j, k, oomph::DenseMatrix< T >::ncol(), oomph::OomphCommunicator::nproc(), oomph::DenseMatrix< T >::nrow(), oomph::LinearAlgebraDistribution::nrow_local(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and p.

Referenced by concatenate_without_communication(), oomph::BlockPreconditioner< MATRIX >::get_block(), oomph::BlockPreconditioner< MATRIX >::get_concatenated_block(), and main().

create_uniformly_distributed_matrix()

void oomph::CRDoubleMatrixHelpers::create_uniformly_distributed_matrix	(	const unsigned &	nrow,
		const unsigned &	ncol,
		const OomphCommunicator *const	comm_pt,
		const Vector< double > &	values,
		const Vector< int > &	column_indices,
		const Vector< int > &	row_start,
		CRDoubleMatrix &	matrix_out
	)

Builds a uniformly distributed matrix. A locally replicated matrix is constructed then redistributed using OOMPH-LIB's default uniform row distribution. This is memory intensive thus should be used for testing or small problems only.

Builds a uniformly distributed matrix. A locally replicated matrix is constructed then redistributed using OOMPH-LIB's default uniform row distribution. This is memory intensive thus should be used for testing or small problems only. The resulting matrix (mat_out) must not have been built.

    {
#ifdef PARANOID
      // Check if the communicator exists.
      if (comm_pt == 0)
      {
        std::ostringstream error_message;
        error_message << "Please supply the communicator.\n";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
      // Is the out matrix built? We need an empty matrix!
      if (matrix_out.built())
      {
        std::ostringstream error_message;
        error_message << "The result matrix has been built.\n"
                      << "Please clear the matrix.\n";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Create the locally replicated distribution.
      bool distributed = false;
      LinearAlgebraDistribution locally_replicated_distribution(
        comm_pt, nrow, distributed);
 
      // Create the matrix.
      matrix_out.build(&locally_replicated_distribution,
                       ncol,
                       values,
                       column_indices,
                       row_start);
 
      // Create the distributed distribution.
      distributed = true;
      LinearAlgebraDistribution distributed_distribution(
        comm_pt, nrow, distributed);
 
      // Redistribute the matrix.
      matrix_out.redistribute(&distributed_distribution);
    }

References oomph::CRDoubleMatrix::build(), oomph::CRDoubleMatrix::built(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and oomph::CRDoubleMatrix::redistribute().

Referenced by main().

deep_copy()

void oomph::CRDoubleMatrixHelpers::deep_copy ( const CRDoubleMatrix *const  in_matrix_pt, CRDoubleMatrix &  out_matrix  )

inline

Create a deep copy of the matrix pointed to by in_matrix_pt.

    {
#ifdef PARANOID
      // Is the out matrix built? We need an empty out matrix!
      if (out_matrix.built())
      {
        std::ostringstream err_msg;
        err_msg << "The result matrix has been built.\n"
                << "Please clear the matrix.\n";
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
 
      // Check that the in matrix pointer is not null.
      if (in_matrix_pt == 0)
      {
        std::ostringstream err_msg;
        err_msg << "The in_matrix_pt is null.\n";
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
 
      // Check that the in matrix is built.
      if (!in_matrix_pt->built())
      {
        std::ostringstream err_msg;
        err_msg << "The in_matrix_pt is null.\n";
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // First set the matrix matrix multiply methods (for both serial and
      // distributed)
      out_matrix.serial_matrix_matrix_multiply_method() =
        in_matrix_pt->serial_matrix_matrix_multiply_method();
 
      out_matrix.distributed_matrix_matrix_multiply_method() =
        in_matrix_pt->distributed_matrix_matrix_multiply_method();
 
 
      // The local nrow and nnz of the in matrix
      const unsigned in_nrow_local = in_matrix_pt->nrow_local();
      const unsigned long in_nnz = in_matrix_pt->nnz();
 
      // Storage for the values, column indices and row start
      double* out_values = new double[in_nnz];
      int* out_column_indices = new int[in_nnz];
      int* out_row_start = new int[in_nrow_local + 1];
 
      // The data to copy over
      const double* const in_values = in_matrix_pt->value();
      const int* const in_column_indices = in_matrix_pt->column_index();
      const int* const in_row_start = in_matrix_pt->row_start();
 
      // Copy the data
      std::copy(in_values, in_values + in_nnz, out_values);
 
      std::copy(
        in_column_indices, in_column_indices + in_nnz, out_column_indices);
 
      std::copy(
        in_row_start, in_row_start + (in_nrow_local + 1), out_row_start);
 
      // Build the matrix
      out_matrix.build(in_matrix_pt->distribution_pt());
 
      out_matrix.build_without_copy(in_matrix_pt->ncol(),
                                    in_nnz,
                                    out_values,
                                    out_column_indices,
                                    out_row_start);
 
      // The only thing we haven't copied over is the default linear solver
      // pointer, but I cannot figure out how to copy over a solver since
      // I do not know what it is.
    } // EoFunc deep_copy

References oomph::CRDoubleMatrix::build(), oomph::CRDoubleMatrix::build_without_copy(), oomph::CRDoubleMatrix::built(), oomph::CRDoubleMatrix::column_index(), copy(), oomph::CRDoubleMatrix::distributed_matrix_matrix_multiply_method(), oomph::DistributableLinearAlgebraObject::distribution_pt(), oomph::CRDoubleMatrix::ncol(), oomph::CRDoubleMatrix::nnz(), oomph::DistributableLinearAlgebraObject::nrow_local(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::CRDoubleMatrix::row_start(), oomph::CRDoubleMatrix::serial_matrix_matrix_multiply_method(), and oomph::CRDoubleMatrix::value().

Referenced by oomph::DoubleMultiVector::DoubleMultiVector(), oomph::BlockPreconditioner< MATRIX >::get_block(), and oomph::BlockPreconditioner< MATRIX >::get_dof_level_block().

gershgorin_eigenvalue_estimate()

double oomph::CRDoubleMatrixHelpers::gershgorin_eigenvalue_estimate

(

const DenseMatrix< CRDoubleMatrix * > &

matrix_pt

)

Calculates the largest Gershgorin disc whilst preserving the sign. Let A be an n by n matrix, with entries aij. For \( i \in \{ 1,...,n \} \) let \( R_i = \sum_{i\neq j} |a_{ij}| \) be the sum of the absolute values of the non-diagonal entries in the i-th row. Let \( D(a_{ii},R_i) \) be the closed disc centered at \( a_{ii} \) with radius \( R_i \), such a disc is called a Gershgorin disc.

We calculate \( |D(a_{ii},R_i)|_{max} \) and multiply by the sign of the diagonal entry.

The DenseMatrix of CRDoubleMatrices are treated as if they are one large matrix. Therefore the dimensions of the sub matrices has to "make sense", there is a paranoid check for this.

Calculates the largest Gershgorin disc whilst preserving the sign. Let A be an n by n matrix, with entries aij. For \( i \in \{ 1,...,n \} \) let \( R_i = \sum_{i\neq j} |a_{ij}| \) be the sum of the absolute values of the non-diagonal entries in the i-th row. Let \( D(a_{ii},R_i) \) be the closed disc centered at \( a_{ii} \) with radius \( R_i \), such a disc is called a Gershgorin disc.

We calculate \( |D(a_{ii},R_i)|_{max} \)and multiply by the sign of the diagonal entry.

The DenseMatrix of CRDoubleMatrices are treated as if they are one large matrix. Therefore the dimensions of the sub matrices has to "make sense", there is a paranoid check for this.

    {
      // The number of block rows and columns
      const unsigned nblockrow = matrix_pt.nrow();
      const unsigned nblockcol = matrix_pt.ncol();
 
#ifdef PARANOID
      // Check that tehre is at least one matrix.
      if (matrix_pt.nrow() == 0)
      {
        std::ostringstream error_message;
        error_message << "There are no matrices... \n";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
 
 
      // Check that all matrix_pt pointers are not null
      // and the matrices are built.
      for (unsigned block_row_i = 0; block_row_i < nblockrow; block_row_i++)
      {
        for (unsigned block_col_i = 0; block_col_i < nblockcol; block_col_i++)
        {
          if (matrix_pt(block_row_i, block_col_i) == 0)
          {
            std::ostringstream error_message;
            error_message << "The pointer martrix_pt(" << block_row_i << ","
                          << block_col_i << ") is null.\n";
            throw OomphLibError(error_message.str(),
                                OOMPH_CURRENT_FUNCTION,
                                OOMPH_EXCEPTION_LOCATION);
          }
 
          if (!matrix_pt(block_row_i, block_col_i)->built())
          {
            std::ostringstream error_message;
            error_message << "The matrix at martrix_pt(" << block_row_i << ","
                          << block_col_i << ") is not built.\n";
            throw OomphLibError(error_message.str(),
                                OOMPH_CURRENT_FUNCTION,
                                OOMPH_EXCEPTION_LOCATION);
          }
        }
      }
#endif
 
 
#ifdef OOMPH_HAS_MPI
 
      // The communicator pointer from block (0,0)
      // All communicators should be the same, we check this next.
      const OomphCommunicator* const comm_pt =
        matrix_pt(0, 0)->distribution_pt()->communicator_pt();
 
#ifdef PARANOID
 
      // Check that all communicators are the same
      for (unsigned block_row_i = 0; block_row_i < nblockrow; block_row_i++)
      {
        for (unsigned block_col_i = 0; block_col_i < nblockcol; block_col_i++)
        {
          // Communicator for this block matrix.
          const OomphCommunicator current_block_comm =
            *(matrix_pt(block_row_i, block_col_i)
                ->distribution_pt()
                ->communicator_pt());
          if (*comm_pt != current_block_comm)
          {
            std::ostringstream error_message;
            error_message << "The communicator of block martrix_pt("
                          << block_row_i << "," << block_col_i
                          << ") is not the same as block "
                          << "matrix_pt(0,0).\n";
            throw OomphLibError(error_message.str(),
                                OOMPH_CURRENT_FUNCTION,
                                OOMPH_EXCEPTION_LOCATION);
          }
        }
      }
 
      // Check that all distributed boolean are the same (if on more than 1
      // core)
      if (comm_pt->nproc() > 1)
      {
        // Get the distributed boolean from matrix_pt(0,0)
        bool first_distributed = matrix_pt(0, 0)->distributed();
 
        for (unsigned block_row_i = 0; block_row_i < nblockrow; block_row_i++)
        {
          for (unsigned block_col_i = 0; block_col_i < nblockcol; block_col_i++)
          {
            // Is the current block distributed?
            bool current_distributed =
              matrix_pt(block_row_i, block_col_i)->distributed();
 
            if (first_distributed != current_distributed)
            {
              std::ostringstream error_message;
              error_message << "Block matrix_pt(" << block_row_i << ","
                            << block_col_i << ") and block matrix_pt(0,0) "
                            << "have a different distributed boolean.\n";
              throw OomphLibError(error_message.str(),
                                  OOMPH_CURRENT_FUNCTION,
                                  OOMPH_EXCEPTION_LOCATION);
            }
          }
        }
      }
 
      // Check that all sub matrix dimensions "make sense"
      // We need to check that all the matrices in the same row has the same
      // nrow. Then repeat for the columns.
 
      // Check the nrow of each block row.
      for (unsigned block_row_i = 0; block_row_i < nblockrow; block_row_i++)
      {
        // Get the nrow to compare against from the first column.
        const unsigned first_block_nrow = matrix_pt(block_row_i, 0)->nrow();
 
        // Loop through the block columns.
        for (unsigned block_col_i = 1; block_col_i < nblockcol; block_col_i++)
        {
          // If the nrow of this block is not the same as the nrow from the
          // first block in this block row, throw an error.
          const unsigned current_block_nrow =
            matrix_pt(block_row_i, block_col_i)->nrow();
 
          if (first_block_nrow != current_block_nrow)
          {
            std::ostringstream error_message;
            error_message << "First block has nrow = " << current_block_nrow
                          << ". But martrix_pt(" << block_row_i << ","
                          << block_col_i
                          << ") has nrow = " << current_block_nrow << ".\n";
            throw OomphLibError(error_message.str(),
                                OOMPH_CURRENT_FUNCTION,
                                OOMPH_EXCEPTION_LOCATION);
          }
        }
      }
 
      // Check the ncol of each block column.
      for (unsigned block_col_i = 0; block_col_i < nblockcol; block_col_i++)
      {
        // Get the ncol from the first block row to compare against.
        const unsigned first_block_ncol = matrix_pt(0, block_col_i)->ncol();
 
        for (unsigned block_row_i = 1; block_row_i < nblockrow; block_row_i++)
        {
          // Get the ncol for the current block.
          const unsigned current_block_ncol =
            matrix_pt(block_row_i, block_col_i)->ncol();
 
          if (first_block_ncol != current_block_ncol)
          {
            std::ostringstream error_message;
            error_message << "First block has ncol = " << current_block_ncol
                          << ". But martrix_pt(" << block_row_i << ","
                          << block_col_i
                          << ") has ncol = " << current_block_ncol << ".\n";
            throw OomphLibError(error_message.str(),
                                OOMPH_CURRENT_FUNCTION,
                                OOMPH_EXCEPTION_LOCATION);
          }
        }
      }
 
      // Check that the distribution for each block row is the same.
      for (unsigned block_row_i = 0; block_row_i < nblockrow; block_row_i++)
      {
        // The first distribution of this block row.
        const LinearAlgebraDistribution first_dist =
          *(matrix_pt(block_row_i, 0)->distribution_pt());
 
        // Loop through the rest of the block columns.
        for (unsigned block_col_i = 1; block_col_i < nblockcol; block_col_i++)
        {
          // Get the distribution from the current block.
          const LinearAlgebraDistribution current_dist =
            matrix_pt(block_row_i, block_col_i)->distribution_pt();
 
          // Compare the first distribution against the current.
          if (first_dist != current_dist)
          {
            std::ostringstream error_message;
            error_message << "First distribution of block row " << block_row_i
                          << " is different from the distribution from "
                          << "martrix_pt(" << block_row_i << "," << block_col_i
                          << ").\n";
            throw OomphLibError(error_message.str(),
                                OOMPH_CURRENT_FUNCTION,
                                OOMPH_EXCEPTION_LOCATION);
          }
        }
      }
 
#endif
#endif
 
      // Loop thrpugh the block rows, then block columns to
      // compute the local inf norm
      double extreme_disc = 0;
      for (unsigned block_row_i = 0; block_row_i < nblockrow; block_row_i++)
      {
        // Get the number of local rows from the first block.
        unsigned block_nrow_local = matrix_pt(block_row_i, 0)->nrow_local();
 
        // Loop through the block_nrow_local in this block row
        for (unsigned local_row_i = 0; local_row_i < block_nrow_local;
             local_row_i++)
        {
          double abs_sum_of_row = 0;
          // Loop through the block columns
          for (unsigned block_col_i = 0; block_col_i < nblockcol; block_col_i++)
          {
            // Locally cache the pointer to the current block.
            CRDoubleMatrix* block_pt = matrix_pt(block_row_i, block_col_i);
 
            const int* row_start = block_pt->row_start();
            const double* value = block_pt->value();
 
            // Loop through the values
            for (int val_i = row_start[local_row_i];
                 val_i < row_start[local_row_i + 1];
                 val_i++)
            {
              abs_sum_of_row += fabs(value[val_i]);
            }
          }
 
          // Now minus the diagonal entry...
          // Locate the diagonal block matrix.
          double* s_values = matrix_pt(block_row_i, block_row_i)->value();
          int* s_column_index =
            matrix_pt(block_row_i, block_row_i)->column_index();
          int* s_row_start = matrix_pt(block_row_i, block_row_i)->row_start();
          // int s_nrow_local =
          // matrix_pt(block_row_i,block_row_i)->nrow_local();
          int s_first_row = matrix_pt(block_row_i, block_row_i)->first_row();
 
          // Get the diagonal value...
          double diagonal_value = 0;
          bool found = false;
          for (int j = s_row_start[local_row_i];
               j < s_row_start[local_row_i + 1] && !found;
               j++)
          {
            if (s_column_index[j] == int(local_row_i + s_first_row))
            {
              diagonal_value = s_values[j];
              found = true;
            }
          }
 
          // Check if the diagonal entry is found.
          if (!found)
          {
            std::ostringstream error_message;
            error_message << "The diagonal entry for the block(" << block_row_i
                          << "," << block_row_i << ")\n"
                          << "on local row " << local_row_i
                          << " does not exist." << std::endl;
            throw OomphLibError(error_message.str(),
                                OOMPH_CURRENT_FUNCTION,
                                OOMPH_EXCEPTION_LOCATION);
          }
 
          // This is the disc.
          abs_sum_of_row -= fabs(diagonal_value);
 
          // Now we have to check if the diagonal entry is
          // on the left or right side of zero.
          if (diagonal_value > 0)
          {
            double extreme_disc_local = diagonal_value + abs_sum_of_row;
            extreme_disc = std::max(extreme_disc_local, extreme_disc);
          }
          else
          {
            double extreme_disc_local = diagonal_value - abs_sum_of_row;
            extreme_disc = std::min(extreme_disc_local, extreme_disc);
          }
        } // Loop through local row (of all block column)
      } // Loop through block row
 
      // if this vector is distributed and on multiple processors then gather
#ifdef OOMPH_HAS_MPI
      double extreme_disc_local = extreme_disc;
      if (matrix_pt(0, 0)->distributed() && comm_pt->nproc() > 1)
      {
        if (extreme_disc > 0)
        {
          MPI_Allreduce(&extreme_disc,
                        &extreme_disc_local,
                        1,
                        MPI_DOUBLE,
                        MPI_MAX,
                        comm_pt->mpi_comm());
        }
        else
        {
          MPI_Allreduce(&extreme_disc,
                        &extreme_disc_local,
                        1,
                        MPI_DOUBLE,
                        MPI_MIN,
                        comm_pt->mpi_comm());
        }
      }
      extreme_disc = extreme_disc_local;
#endif
 
      // and return
      return extreme_disc;
    }

References boost::multiprecision::fabs(), MergeRestartFiles::found, j, max, min, oomph::DenseMatrix< T >::ncol(), oomph::OomphCommunicator::nproc(), oomph::DenseMatrix< T >::nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::CRDoubleMatrix::row_start(), Eigen::value, and oomph::CRDoubleMatrix::value().

inf_norm()

double oomph::CRDoubleMatrixHelpers::inf_norm

(

const DenseMatrix< CRDoubleMatrix * > &

matrix_pt

)

Compute infinity (maximum) norm of sub blocks as if it was one matrix.

Calculates the infinity (maximum) norm of a DenseMartrix of CRDoubleMatrices as if it was one large matrix. This avoids creating a concatenation of the sub-blocks just to calculate the infinity norm.

    {
      // The number of block rows and columns
      const unsigned nblockrow = matrix_pt.nrow();
      const unsigned nblockcol = matrix_pt.ncol();
 
#ifdef PARANOID
      // Check that tehre is at least one matrix.
      if (matrix_pt.nrow() == 0)
      {
        std::ostringstream error_message;
        error_message << "There are no matrices... \n";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
 
 
      // Check that all matrix_pt pointers are not null
      // and the matrices are built.
      for (unsigned block_row_i = 0; block_row_i < nblockrow; block_row_i++)
      {
        for (unsigned block_col_i = 0; block_col_i < nblockcol; block_col_i++)
        {
          if (matrix_pt(block_row_i, block_col_i) == 0)
          {
            std::ostringstream error_message;
            error_message << "The pointer martrix_pt(" << block_row_i << ","
                          << block_col_i << ") is null.\n";
            throw OomphLibError(error_message.str(),
                                OOMPH_CURRENT_FUNCTION,
                                OOMPH_EXCEPTION_LOCATION);
          }
 
          if (!matrix_pt(block_row_i, block_col_i)->built())
          {
            std::ostringstream error_message;
            error_message << "The matrix at martrix_pt(" << block_row_i << ","
                          << block_col_i << ") is not built.\n";
            throw OomphLibError(error_message.str(),
                                OOMPH_CURRENT_FUNCTION,
                                OOMPH_EXCEPTION_LOCATION);
          }
        }
      }
#endif
 
#ifdef OOMPH_HAS_MPI
 
      // The communicator pointer from block (0,0)
      const OomphCommunicator* const comm_pt =
        matrix_pt(0, 0)->distribution_pt()->communicator_pt();
 
#ifdef PARANOID
 
 
      // Check that all communicators are the same
      for (unsigned block_row_i = 0; block_row_i < nblockrow; block_row_i++)
      {
        for (unsigned block_col_i = 0; block_col_i < nblockcol; block_col_i++)
        {
          // Communicator for this block matrix.
          const OomphCommunicator current_block_comm =
            *(matrix_pt(block_row_i, block_col_i)
                ->distribution_pt()
                ->communicator_pt());
          if (*comm_pt != current_block_comm)
          {
            std::ostringstream error_message;
            error_message << "The communicator of block martrix_pt("
                          << block_row_i << "," << block_col_i
                          << ") is not the same as block "
                          << "matrix_pt(0,0).\n";
            throw OomphLibError(error_message.str(),
                                OOMPH_CURRENT_FUNCTION,
                                OOMPH_EXCEPTION_LOCATION);
          }
        }
      }
 
      // Check that all distributed boolean are the same (if on more than 1
      // core)
      if (comm_pt->nproc() > 1)
      {
        // Get the distributed boolean from matrix_pt(0,0)
        bool first_distributed = matrix_pt(0, 0)->distributed();
 
        for (unsigned block_row_i = 0; block_row_i < nblockrow; block_row_i++)
        {
          for (unsigned block_col_i = 0; block_col_i < nblockcol; block_col_i++)
          {
            // Is the current block distributed?
            bool current_distributed =
              matrix_pt(block_row_i, block_col_i)->distributed();
 
            if (first_distributed != current_distributed)
            {
              std::ostringstream error_message;
              error_message << "Block matrix_pt(" << block_row_i << ","
                            << block_col_i << ") and block matrix_pt(0,0) "
                            << "have a different distributed boolean.\n";
              throw OomphLibError(error_message.str(),
                                  OOMPH_CURRENT_FUNCTION,
                                  OOMPH_EXCEPTION_LOCATION);
            }
          }
        }
      }
 
      // Check that all sub matrix dimensions "make sense"
      // We need to check that all the matrices in the same row has the same
      // nrow. Then repeat for the columns.
 
      // Check the nrow of each block row.
      for (unsigned block_row_i = 0; block_row_i < nblockrow; block_row_i++)
      {
        // Get the nrow to compare against from the first column.
        const unsigned first_block_nrow = matrix_pt(block_row_i, 0)->nrow();
 
        // Loop through the block columns.
        for (unsigned block_col_i = 1; block_col_i < nblockcol; block_col_i++)
        {
          // If the nrow of this block is not the same as the nrow from the
          // first block in this block row, throw an error.
          const unsigned current_block_nrow =
            matrix_pt(block_row_i, block_col_i)->nrow();
 
          if (first_block_nrow != current_block_nrow)
          {
            std::ostringstream error_message;
            error_message << "First block has nrow = " << current_block_nrow
                          << ". But martrix_pt(" << block_row_i << ","
                          << block_col_i
                          << ") has nrow = " << current_block_nrow << ".\n";
            throw OomphLibError(error_message.str(),
                                OOMPH_CURRENT_FUNCTION,
                                OOMPH_EXCEPTION_LOCATION);
          }
        }
      }
 
      // Check the ncol of each block column.
      for (unsigned block_col_i = 0; block_col_i < nblockcol; block_col_i++)
      {
        // Get the ncol from the first block row to compare against.
        const unsigned first_block_ncol = matrix_pt(0, block_col_i)->ncol();
 
        for (unsigned block_row_i = 1; block_row_i < nblockrow; block_row_i++)
        {
          // Get the ncol for the current block.
          const unsigned current_block_ncol =
            matrix_pt(block_row_i, block_col_i)->ncol();
 
          if (first_block_ncol != current_block_ncol)
          {
            std::ostringstream error_message;
            error_message << "First block has ncol = " << current_block_ncol
                          << ". But martrix_pt(" << block_row_i << ","
                          << block_col_i
                          << ") has ncol = " << current_block_ncol << ".\n";
            throw OomphLibError(error_message.str(),
                                OOMPH_CURRENT_FUNCTION,
                                OOMPH_EXCEPTION_LOCATION);
          }
        }
      }
 
      // Check that the distribution for each block row is the same.
      for (unsigned block_row_i = 0; block_row_i < nblockrow; block_row_i++)
      {
        // The first distribution of this block row.
        const LinearAlgebraDistribution first_dist =
          *(matrix_pt(block_row_i, 0)->distribution_pt());
 
        // Loop through the rest of the block columns.
        for (unsigned block_col_i = 1; block_col_i < nblockcol; block_col_i++)
        {
          // Get the distribution from the current block.
          const LinearAlgebraDistribution current_dist =
            matrix_pt(block_row_i, block_col_i)->distribution_pt();
 
          // Compare the first distribution against the current.
          if (first_dist != current_dist)
          {
            std::ostringstream error_message;
            error_message << "First distribution of block row " << block_row_i
                          << " is different from the distribution from "
                          << "martrix_pt(" << block_row_i << "," << block_col_i
                          << ").\n";
            throw OomphLibError(error_message.str(),
                                OOMPH_CURRENT_FUNCTION,
                                OOMPH_EXCEPTION_LOCATION);
          }
        }
      }
#endif
 
#endif
 
      // Loop thrpugh the block rows, then block columns to
      // compute the local inf norm
      double inf_norm = 0;
      for (unsigned block_row_i = 0; block_row_i < nblockrow; block_row_i++)
      {
        // Get the number of local rows from the first block.
        unsigned block_nrow_local = matrix_pt(block_row_i, 0)->nrow_local();
 
        // Loop through the block_nrow_local in this block row
        for (unsigned local_row_i = 0; local_row_i < block_nrow_local;
             local_row_i++)
        {
          double abs_sum_of_row = 0;
          // Loop through the block columns
          for (unsigned block_col_i = 0; block_col_i < nblockcol; block_col_i++)
          {
            // Locally cache the pointer to the current block.
            CRDoubleMatrix* block_pt = matrix_pt(block_row_i, block_col_i);
 
            const int* row_start = block_pt->row_start();
            const double* value = block_pt->value();
 
            // Loop through the values
            for (int val_i = row_start[local_row_i];
                 val_i < row_start[local_row_i + 1];
                 val_i++)
            {
              abs_sum_of_row += fabs(value[val_i]);
            }
          }
          // Store the max row
          inf_norm = std::max(inf_norm, abs_sum_of_row);
        }
      }
 
      // if this vector is distributed and on multiple processors then gather
#ifdef OOMPH_HAS_MPI
      double inf_norm_local = inf_norm;
      if (matrix_pt(0, 0)->distributed() && comm_pt->nproc() > 1)
      {
        MPI_Allreduce(&inf_norm,
                      &inf_norm_local,
                      1,
                      MPI_DOUBLE,
                      MPI_MAX,
                      comm_pt->mpi_comm());
      }
      inf_norm = inf_norm_local;
#endif
 
      // and return
      return inf_norm;
    }

References boost::multiprecision::fabs(), max, oomph::DenseMatrix< T >::ncol(), oomph::OomphCommunicator::nproc(), oomph::DenseMatrix< T >::nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::CRDoubleMatrix::row_start(), Eigen::value, and oomph::CRDoubleMatrix::value().

Referenced by oomph::PseudoElasticPreconditioner::setup(), and oomph::LagrangeEnforcedFlowPreconditioner::setup().