Отчёт_неявная_схема (Весенний семестр), страница 6

Westill keep it, since other processors will be busy with them anyway.alpha_coeff_x_axis[0] = 0;if (coords[0] == 0) {F_coeff_x_axis[0] -= A_coeff_x_axis * upper_border[i];alpha_coeff_x_axis[0] = 0;up_data_buffer[0] = 0;}if (coords[0] == dim_sizes[0] - 1) {F_coeff_x_axis[field_length - 1] -= B_coeff_x_axis *34lower_border[i];beta_coeff_x_axis[field_length - 1] = 0;}// Now we have to solve the equation system, which is done in severalbasic steps.//------------------------------------------------------------------------------// At first step, we subtract equation lines top-down; this process makesall alpha coeffs = 0 and adjusts f_coeffs.

However,// It also causes up data buffers to become non-zero and thus makingthose "defective" lines.for (int j = 1; j < field_length; ++j) {// Subtraction coefficient calculation.subtract_coeff_x_axis = alpha_coeff_x_axis[j] /charlie_coeff_x_axis[j - 1];// Solution vector.F_coeff_x_axis[j] -= subtract_coeff_x_axis * F_coeff_x_axis[j 1];// Alpha and charlie vectors.

Beta vector doesn't change duringthis step.alpha_coeff_x_axis[j] = 0; // -= subtract_coeff_x_axis *charlie_coeff_x_axis[j - 1];charlie_coeff_x_axis[j] -= beta_coeff_x_axis[j - 1] *subtract_coeff_x_axis;// Up data buffer gets filled with defective data :(up_data_buffer[j] -= up_data_buffer[j - 1] *subtract_coeff_x_axis;}// At second step same sequence of operations is done in reversal.// This nullifies beta coefficient and fills down data buffers withunnecessary data.// Down data buffer values are initialized here;if (field_length >= 2) {down_data_buffer[field_length - 1] =charlie_coeff_x_axis[field_length - 1];down_data_buffer[field_length - 2] =beta_coeff_x_axis[field_length - 2];}// This also changes data in up data buffer !!!for (int j = field_length - 3; j >= 0; --j) {// Subtraction coefficient calculation.subtract_coeff_x_axis = beta_coeff_x_axis[j] /charlie_coeff_x_axis[j + 1];// Solution vector.F_coeff_x_axis[j] -= subtract_coeff_x_axis * F_coeff_x_axis[j +1];35// Beta vector becomes zero, charlie vector is not affected( alpha is zero ).beta_coeff_x_axis[j] -= subtract_coeff_x_axis *charlie_coeff_x_axis[j + 1];// Down data buffer is filled with defective data.

Up databuffer is adjusted.up_data_buffer[j] -= up_data_buffer[j + 1] *subtract_coeff_x_axis;down_data_buffer[j] -= down_data_buffer[j + 1] *subtract_coeff_x_axis;}// Now the tricky part: we take the upper line, take CHARLIE, UP, DOWNand SOLUTION// values from there and send them to the previous process, if there'sone.send_buf[0]send_buf[1]send_buf[2]send_buf[3]====up_data_buffer[0];charlie_coeff_x_axis[0];down_data_buffer[0];F_coeff_x_axis[0];if (dim_sizes[0] > 1) { // There is somebody to receive.MPI_Cart_shift(cart_comm, 0, -1, &shift_source, &shift_dest);MPI_Sendrecv(send_buf, 4, MPI_DOUBLE, shift_dest, 0, recv_buf,4, MPI_DOUBLE, shift_source, 0, cart_comm, &status);}else {for (int k = 0; k <= 3; ++k) { recv_buf[k] = send_buf[k]; }}// Then we adjust actual values according to what we received and fillONE line of a global matrix.if (coords[0] != dim_sizes[0] - 1) { // It's not the "lowest" process.// Calculate subtraction coefficient.subtract_coeff_x_axis = beta_coeff_x_axis[field_length - 1] /recv_buf[1];// Calculate actual new values.charlie_coeff_x_axis[field_length - 1] -= subtract_coeff_x_axis* recv_buf[0];beta_coeff_x_axis[field_length - 1] = 0;appendix -= subtract_coeff_x_axis * recv_buf[2];F_coeff_x_axis[field_length - 1] -= subtract_coeff_x_axis *recv_buf[3];}// Form a line for the global matrix.send_buf[0] = up_data_buffer[field_length - 1];send_buf[1] = charlie_coeff_x_axis[field_length - 1];send_buf[2] = appendix;send_buf[3] = F_coeff_x_axis[field_length - 1];36if (dim_sizes[0] > 1) { // Send only there is somebody to receive.if (coords[0] == 0) { // Now processes with coords[0] == 0receive data to global matrix.dub_coords[1] = coords[1];for (int j = 0; j <= 3; ++j) { global_x_mtrx[j] =send_buf[j]; }for (int j = 1; j < dim_sizes[0]; ++j) {dub_coords[0] = j;MPI_Cart_rank(cart_comm, dub_coords,&shift_source);MPI_Recv(&global_x_mtrx[4 * j], 4, MPI_DOUBLE,shift_source, 0, cart_comm, &status);}}else { // And other processes send data to them.dub_coords[0] = 0;dub_coords[1] = coords[1];MPI_Cart_rank(cart_comm, dub_coords, &shift_dest);MPI_Send(send_buf, 4, MPI_DOUBLE, shift_dest, 0,cart_comm);}}else {for (int k = 0; k <= 3; ++k) { global_x_mtrx[k] = send_buf[k]; }}// Simple version of the method is used for solving the equation withglobal matrix.simple_solve(global_x_mtrx, global_x_array, dim_sizes[0]);// Results are sent to processes, directly into appopriate locations oftheir dub_arrays.if (dim_sizes[0] > 1) { // There is somebody to receive.if (coords[0] == 0) {array_data[(field_width - 1) * field_length + i] =global_x_array[0];dub_coords[1] = coords[1];for (int j = 1; j < dim_sizes[1]; ++j) {dub_coords[0] = j;MPI_Cart_rank(cart_comm, dub_coords,&shift_source);MPI_Send(&global_x_array[j], 1, MPI_DOUBLE,shift_source, 2, cart_comm);MPI_Send(&global_x_array[j - 1], 1, MPI_DOUBLE,shift_source, 2, cart_comm);}}else {MPI_Recv(&array_data[(field_width - 1) * field_length +i], 1, MPI_DOUBLE, shift_dest,2, cart_comm, &status);MPI_Recv(&auxillary, 1, MPI_DOUBLE, shift_dest,2, cart_comm, &status);37}}else {array_data[(field_width - 1) * field_length + i] =global_x_array[0];}// Finally, every process locally fills the remaining dub_array cells.for (int j = field_length - 2; j >= 0; --j) {array_data[j * field_length + i] =(F_coeff_x_axis[j] - up_data_buffer[j] * auxillary array_data[(field_width - 1) * field_length + i] * down_data_buffer[j]) /charlie_coeff_x_axis[j];}}// As horizontal solution is over, we should have fully updated data ready forthe next step.

That means iteration is over.// --------------------------------------------------------------// --------------------------------------------------------------}// Output results to a file.void Trans_solver::Output(const int mode) {double send_double = 0;double recv_double = 0;// Calculate source process.if (!((coords[0] == 0) && (coords[1] == 0))) { // Need to receive.dub_coords[0] = coords[0];dub_coords[1] = coords[1] - 1;if (dub_coords[1] < 0) {dub_coords[1] = dim_sizes[1] - 1;dub_coords[0]--;}MPI_Cart_rank(cart_comm, dub_coords, &shift_source);MPI_Recv(&recv_double, 1, MPI_DOUBLE, shift_source, 5, cart_comm,&status);}std::fstream out_file;// In parallel version processes write into memory one after another.if (mode == 1) { // Use stdout.for (int i = 0; i < field_width; ++i) {for (int j = 0; j < field_length; ++j) {std::cout << array_data[i * field_length + j] << " ";}}}// In parallel version processes write into memory one after another.else { // Use fstream.38if ((coords[0] == 0) && (coords[1] == 0)) {out_file.open("std_out.txt", std::fstream::out);}else {out_file.open("std_out.txt", std::fstream::in |std::fstream::out | std::fstream::ate);}for (int i = 0; i < field_length; ++i) {for (int j = 0; j < field_width; ++j) {out_file << array_data[i * field_width + j] << " ";}}out_file.close();}// Calculate dest_process.if (!((coords[0] == dim_sizes[0] - 1) && (coords[1] == dim_sizes[1] - 1))) { //Need to send.dub_coords[1] = coords[1] + 1;dub_coords[0] = coords[0];if (dub_coords[1] >= dim_sizes[1]) {dub_coords[1] = 0;dub_coords[0]++;}MPI_Cart_rank(cart_comm, dub_coords, &shift_dest);MPI_Send(&send_double, 1, MPI_DOUBLE, shift_dest, 5, cart_comm);}}// Default destructor.Trans_solver::~Trans_solver(){// Delete cart data.if (dim_sizes != NULL) { delete dim_sizes; }if (periods != NULL) { delete periods; }if (coords != NULL) { delete coords; }if (dub_coords != NULL) { delete dub_coords; }// Delete field data.if (array_data != NULL) { delete array_data; }if (dub_array != NULL) { delete dub_array; }// Delete auxillary data fields.if (left_data_buffer != NULL) { delete left_data_buffer; }if (right_data_buffer != NULL) { delete right_data_buffer; }//ififififDelete border data.(upper_border != NULL) { delete upper_border; }(lower_border != NULL) { delete lower_border; }(left_border != NULL) { delete left_border; }(right_border != NULL) { delete right_border; }// Delete tridiagonal coeffs data.if (alpha_coeff_x_axis != NULL) { delete alpha_coeff_x_axis; }39if (beta_coeff_x_axis != NULL) { delete beta_coeff_x_axis; }if (charlie_coeff_x_axis != NULL) { delete charlie_coeff_x_axis; }if (alpha_coeff_y_axis != NULL) { delete alpha_coeff_y_axis; }if (beta_coeff_y_axis != NULL) { delete beta_coeff_y_axis; }if (charlie_coeff_y_axis != NULL) { delete charlie_coeff_y_axis; }// Delete coefficients data.if (F_coeff_x_axis != NULL) { delete F_coeff_x_axis; }if (F_coeff_y_axis != NULL) { delete F_coeff_y_axis; }// Delete sendrecv buffers data.if (send_buf != NULL) { delete send_buf; }if (recv_buf != NULL) { delete recv_buf; }// Delete global matrices data.if (global_x_mtrx != NULL) { delete global_x_mtrx; }if (global_y_mtrx != NULL) { delete global_y_mtrx; }// Delete global arrays data.if (global_x_array != NULL) { delete global_x_array; }if (global_y_array != NULL) { delete global_y_array; }}// Basic#include#include#includeinclusions.<cstdlib><cstdio><iostream>// Mpi inclusions.#include "mpi.h"// Using directives.using namespace std;// Project inclusions.#include "Trans_solver.cpp"// Main function :)// ------------------------------------------------------------------------------int main(int argc, char* argv[]) {// MPI rank & size variables.int p_rank, p_size;// Command line arguments.int data_length = atoi(argv[1]);int data_width = atoi(argv[2]);int iter_num = atoi(argv[3]);int calc_mode = atoi(argv[4]);int out_mode = atoi(argv[5]);// Initialize MPI.MPI_Init(&argc, &argv);MPI_Comm_rank(MPI_COMM_WORLD, &p_rank);40MPI_Comm_size(MPI_COMM_WORLD, &p_size);// Time parameters.double start_time;double end_time;Trans_solver Eq_solver; // Create a solver.Eq_solver.Init(data_length, data_width, p_rank, p_size); // Initialize solver.Eq_solver.Fill_data(0); // Load test data to solver.start_time = MPI_Wtime();// Perform as many iterations as needed.// Perform as many iterations as needed.if (calc_mode == 0) { // Calculation mode.for (int i = 0; i < iter_num; ++i) {Eq_solver.Iterate(i);}}else if (calc_mode == 1) { // Generate test solutionEq_solver.Fill_data(iter_num);}MPI_Barrier(MPI_COMM_WORLD);end_time = MPI_Wtime();double total_time = end_time - start_time;if (out_mode == 0) { // Calculation mode.if (p_rank == 0) {cout << "Calculation time is : " << total_time << "\n";}}else { // VerifyEq_solver.Output(out_mode);}// Finalize MPI.MPI_Finalize();// Return :)return 0;}41.