diff --git a/raptor/util/linalg/relax.cpp b/raptor/util/linalg/relax.cpp index ad412fe2..ed8d40b1 100644 --- a/raptor/util/linalg/relax.cpp +++ b/raptor/util/linalg/relax.cpp @@ -138,21 +138,23 @@ void ssor(CSRMatrix* A, Vector& b, Vector& x, Vector& tmp, int num_sweeps, // Inverts block at address A // Returns inverted block at address A_inv -void invert_block(double* A, double* D_inv, int n) +void invert_block(double* A, int n) { int info; - int *lu = new int[n]; int block_size = n*n; + int *lu = new int[n]; + double* tmp = new double[block_size]; dgetrf_(&n, &n, A, &n, lu, &info); - dgetri_(&n, A, &n, lu, D_inv, &block_size, &info); + dgetri_(&n, A, &n, lu, tmp, &block_size, &info); + delete[] tmp; delete[] lu; } void block_relax_init(BSRMatrix* A, double** D_inv_ptr) { - double* D_inv = new double[A->n_rows*A->b_size]; + double* D_inv = new double[A->n_rows*A->b_size](); double** bdata = (double**)(A->get_data()); int row_start, row_end; @@ -167,7 +169,9 @@ void block_relax_init(BSRMatrix* A, double** D_inv_ptr) int col = A->idx2[j]; if (i == col) { - invert_block(bdata[j], &(D_inv[i*A->b_size]), A->b_rows); + memcpy(&(D_inv[i*A->b_size]), bdata[j], A->b_size*sizeof(double)); + invert_block(&(D_inv[i*A->b_size]), A->b_rows); + break; } } @@ -181,12 +185,12 @@ void block_relax_init(BSRMatrix* A, double** D_inv_ptr) void jacobi(BSRMatrix* A, double* D_inv, Vector& b, Vector& x, Vector& tmp, int num_sweeps, double omega) { - double* rsum = new double[A->b_size]; - double* tmp_rsum = new double[A->b_size]; + double* rsum = new double[A->b_rows]; + double* tmp_rsum = new double[A->b_rows]; double** bdata = (double**)(A->get_data()); int row_start, row_end; - char trans = 'N'; + char trans = 'T'; int n = A->b_rows; int incr = 1; double one = 1.0; @@ -196,7 +200,7 @@ void jacobi(BSRMatrix* A, double* D_inv, Vector& b, Vector& x, Vector& tmp, int for (int iter = 0; iter < num_sweeps; iter++) { // Copy x to tmp vector - memcpy(tmp.data(), x.data(), tmp.size()); + memcpy(tmp.data(), x.data(), tmp.size()*sizeof(double)); // Begin block Jacobi sweep for (int row = 0; row < A->n_rows; row++) @@ -205,9 +209,7 @@ void jacobi(BSRMatrix* A, double* D_inv, Vector& b, Vector& x, Vector& tmp, int row_end = A->idx1[row+1]; if (row_start == row_end) continue; - int b_row_idx = row * A->b_rows; - - memset(rsum, 0, A->b_size*sizeof(double)); + memset(rsum, 0, n*sizeof(double)); // Block dot product between block row and vector x for (int j = row_start; j < row_end; j++) @@ -215,23 +217,25 @@ void jacobi(BSRMatrix* A, double* D_inv, Vector& b, Vector& x, Vector& tmp, int int col = A->idx2[j]; if (row != col) //ignore diagonal { - dgemv_(&trans, &n, &n, &one, bdata[j], &n, &((tmp[col * A->b_cols])), &incr, &zero, tmp_rsum, &incr); + dgemv_(&trans, &n, &n, &one, bdata[j], &n, &((tmp[col * n])), &incr, &one, rsum, &incr); - for (int k = 0; k < A->b_rows; k++) - rsum[k] += tmp_rsum[k]; + //for (int k = 0; k < n; k++) + // rsum[k] += tmp_rsum[k]; } } // r = b - r / diag // in block form, calculate as: block_r = (b - block_r)*D_inv // - for (int k = 0; k < A->b_rows; k++) - rsum[k] = b[b_row_idx + k] - rsum[k]; - dgemv_(&trans, &n, &n, &one, &(D_inv[row*A->b_size]), &n, &(rsum[0]), &incr, &zero, tmp_rsum, &incr); + for (int k = 0; k < n; k++) + rsum[k] = b[row*n + k] - rsum[k]; + dgemv_(&trans, &n, &n, &one, &(D_inv[row*A->b_size]), &n, rsum, &incr, &zero, tmp_rsum, &incr); // Weighted Jacobi calculation for row for (int k = 0; k < A->b_rows; k++) - x[b_row_idx + k] = omega*tmp_rsum[k] + (1.0-omega)*tmp[b_row_idx + k]; + x[row*n + k] = omega*tmp_rsum[k] + (1.0-omega)*tmp[row*n + k]; + //for (int k = 0; k < A->b_size; k++) + // printf("Dinv[%d][%d] = %e\n", row, k, D_inv[row*A->b_size+k]); } } @@ -292,6 +296,7 @@ void sor(BSRMatrix* A, double* D_inv, Vector& b, Vector& x, Vector& tmp, for (int k = 0; k < A->b_rows; k++) x[b_row_idx + k] = omega*tmp_rsum[k] + (1.0-omega)*x[b_row_idx + k]; + } } diff --git a/raptor/util/linalg/relax.hpp b/raptor/util/linalg/relax.hpp index 7bc7995b..6f87fc37 100644 --- a/raptor/util/linalg/relax.hpp +++ b/raptor/util/linalg/relax.hpp @@ -20,14 +20,14 @@ void ssor(CSRMatrix* A, Vector& b, Vector& x, Vector& tmp, int num_sweeps = 1, double omega = 1.0); // Block Methods (BSR Matrix) -void relax_init(BSRMatrix* A, double** D_inv_ptr, int n); +void block_relax_init(BSRMatrix* A, double** D_inv_ptr); void jacobi(BSRMatrix* A, double* D_inv, Vector& b, Vector& x, Vector& tmp, int num_sweeps = 1, double omega = 1.0); void sor(BSRMatrix* A, double* D_inv, Vector& b, Vector& x, Vector& tmp, int num_sweeps = 1, double omega = 1.0); void ssor(BSRMatrix* A, double* D_inv, Vector& b, Vector& x, Vector& tmp, int num_sweeps = 1, double omega = 1.0); -void relax_free(double* D_inv); +void block_relax_free(double* D_inv); } diff --git a/raptor/util/tests/CMakeLists.txt b/raptor/util/tests/CMakeLists.txt index 7633dd88..4266b35c 100644 --- a/raptor/util/tests/CMakeLists.txt +++ b/raptor/util/tests/CMakeLists.txt @@ -54,6 +54,14 @@ add_test(LaplaceGSTest ./test_gs_laplacian) #target_link_libraries(test_bsr_jacobi_aniso raptor ${MPI_LIBRARIES} googletest pthread ) #add_test(BSRAnisoJacobiTest ./test_bsr_jacobi_aniso) +add_executable(test_bsr_jacobi_aniso test_bsr_jacobi_aniso.cpp) +target_link_libraries(test_bsr_jacobi_aniso raptor ${MPI_LIBRARIES} googletest pthread ) +add_test(BSRAnisoJacobiTest ./test_bsr_jacobi_aniso) + +add_executable(test_bsr_jacobi_laplacian test_bsr_jacobi_laplacian.cpp) +target_link_libraries(test_bsr_jacobi_laplacian raptor ${MPI_LIBRARIES} googletest pthread ) +add_test(BSRLaplaceJacobiTest ./test_bsr_jacobi_laplacian) + if (WITH_MPI) add_executable(test_par_add test_par_add.cpp) diff --git a/raptor/util/tests/test_bsr_jacobi_aniso.cpp b/raptor/util/tests/test_bsr_jacobi_aniso.cpp index 99aa884d..21a6eda5 100644 --- a/raptor/util/tests/test_bsr_jacobi_aniso.cpp +++ b/raptor/util/tests/test_bsr_jacobi_aniso.cpp @@ -29,10 +29,16 @@ TEST(AnisoJacobiTest, TestsInUtil) Vector b(A_sten->n_rows); Vector tmp(A_sten->n_rows); - const char* x_ones_1 = "../../../../test_data/aniso_jacobi_ones_1.txt"; - const char* x_ones_2 = "../../../../test_data/aniso_jacobi_ones_2.txt"; - const char* x_inc_1 = "../../../../test_data/aniso_jacobi_inc_1.txt"; - const char* x_inc_2 = "../../../../test_data/aniso_jacobi_inc_2.txt"; + BSRMatrix* A = new BSRMatrix(A_sten, 5, 5); + double* D_inv; + + const char* x_ones_1 = "../../../../test_data/aniso_block_jacobi_ones_1.txt"; + const char* x_ones_2 = "../../../../test_data/aniso_block_jacobi_ones_2.txt"; + const char* x_inc_1 = "../../../../test_data/aniso_block_jacobi_inc_1.txt"; + const char* x_inc_2 = "../../../../test_data/aniso_block_jacobi_inc_2.txt"; + + // SETUP D_inv! + block_relax_init(A, &D_inv); /********************************************* * Test jacobi when b is constant value 1.0 @@ -40,7 +46,7 @@ TEST(AnisoJacobiTest, TestsInUtil) // Iteration 1 b.set_const_value(1.0); x.set_const_value(0.0); - jacobi(A_sten, b, x, tmp, 1, 3.0/4); + jacobi(A, D_inv, b, x, tmp, 1, 3.0/4); f = fopen(x_ones_1, "r"); for (int i = 0; i < A_sten->n_rows; i++) { @@ -51,7 +57,7 @@ TEST(AnisoJacobiTest, TestsInUtil) fclose(f); // Iteration 2 - jacobi(A_sten, b, x, tmp, 1, 3.0/4); + jacobi(A, D_inv, b, x, tmp, 1, 3.0/4); f = fopen(x_ones_2, "r"); for (int i = 0; i < A_sten->n_rows; i++) { @@ -60,6 +66,7 @@ TEST(AnisoJacobiTest, TestsInUtil) ASSERT_NEAR(x[i],x_val,1e-06); } fclose(f); + /********************************************* * Test jacobi when b[i] = i @@ -70,7 +77,7 @@ TEST(AnisoJacobiTest, TestsInUtil) b[i] = i; } x.set_const_value(0.0); - jacobi(A_sten, b, x, tmp, 1, 3.0/4); + jacobi(A, D_inv, b, x, tmp, 1, 3.0/4); f = fopen(x_inc_1, "r"); for (int i = 0; i < A_sten->n_rows; i++) { @@ -81,7 +88,7 @@ TEST(AnisoJacobiTest, TestsInUtil) fclose(f); // Iteration 2 - jacobi(A_sten, b, x, tmp, 1, 3.0/4); + jacobi(A, D_inv, b, x, tmp, 1, 3.0/4); f = fopen(x_inc_2, "r"); for (int i = 0; i < A_sten->n_rows; i++) { @@ -90,7 +97,13 @@ TEST(AnisoJacobiTest, TestsInUtil) ASSERT_NEAR(x[i],x_val,1e-06); } fclose(f); + + + // Free D_inv + block_relax_free(D_inv); + delete A; + delete A_sten; } // end of TEST(AnisoSpMVTest, TestsInUtil) // diff --git a/raptor/util/tests/test_bsr_jacobi_laplacian.cpp b/raptor/util/tests/test_bsr_jacobi_laplacian.cpp new file mode 100644 index 00000000..e4d76b31 --- /dev/null +++ b/raptor/util/tests/test_bsr_jacobi_laplacian.cpp @@ -0,0 +1,108 @@ +// Copyright (c) 2015-2017, RAPtor Developer Team +// License: Simplified BSD, http://opensource.org/licenses/BSD-2-Clause + + +#include "gtest/gtest.h" +#include "raptor/raptor.hpp" + +using namespace raptor; + +int main(int argc, char** argv) +{ + ::testing::InitGoogleTest(&argc, argv); + return RUN_ALL_TESTS(); + +} // end of main() // + +TEST(AnisoJacobiTest, TestsInUtil) +{ + double x_val; + int grid[3] = {10, 10, 10}; + double* stencil = laplace_stencil_27pt(); + CSRMatrix* A_sten = stencil_grid(stencil, grid, 3); + int n_items_read; + FILE* f; + + Vector x(A_sten->n_rows); + Vector b(A_sten->n_rows); + Vector tmp(A_sten->n_rows); + + BSRMatrix* A = new BSRMatrix(A_sten, 5, 5); + double* D_inv; + + const char* x_ones_1 = "../../../../test_data/laplace_block_jacobi_ones_1.txt"; + const char* x_ones_2 = "../../../../test_data/laplace_block_jacobi_ones_2.txt"; + const char* x_inc_1 = "../../../../test_data/laplace_block_jacobi_inc_1.txt"; + const char* x_inc_2 = "../../../../test_data/laplace_block_jacobi_inc_2.txt"; + + // SETUP D_inv! + block_relax_init(A, &D_inv); + + /********************************************* + * Test jacobi when b is constant value 1.0 + *********************************************/ + // Iteration 1 + b.set_const_value(1.0); + x.set_const_value(0.0); + jacobi(A, D_inv, b, x, tmp, 1, 3.0/4); + f = fopen(x_ones_1, "r"); + for (int i = 0; i < A_sten->n_rows; i++) + { + n_items_read = fscanf(f, "%lg\n", &x_val); + ASSERT_EQ(n_items_read, 1); + ASSERT_NEAR(x[i],x_val,1e-06); + } + fclose(f); + + // Iteration 2 + jacobi(A, D_inv, b, x, tmp, 1, 3.0/4); + f = fopen(x_ones_2, "r"); + for (int i = 0; i < A_sten->n_rows; i++) + { + n_items_read = fscanf(f, "%lg\n", &x_val); + ASSERT_EQ(n_items_read, 1); + ASSERT_NEAR(x[i],x_val,1e-06); + } + fclose(f); + + /********************************************* + * Test jacobi when b[i] = i + *********************************************/ + // Iteration 1 + for (int i = 0; i < A_sten->n_rows; i++) + { + b[i] = i; + } + x.set_const_value(0.0); + jacobi(A, D_inv, b, x, tmp, 1, 3.0/4); + f = fopen(x_inc_1, "r"); + for (int i = 0; i < A_sten->n_rows; i++) + { + n_items_read = fscanf(f, "%lg\n", &x_val); + ASSERT_EQ(n_items_read, 1); + ASSERT_NEAR(x[i],x_val,1e-06); + } + fclose(f); + + // Iteration 2 + jacobi(A, D_inv, b, x, tmp, 1, 3.0/4); + f = fopen(x_inc_2, "r"); + for (int i = 0; i < A_sten->n_rows; i++) + { + n_items_read = fscanf(f, "%lg\n", &x_val); + ASSERT_EQ(n_items_read, 1); + ASSERT_NEAR(x[i],x_val,1e-06); + } + fclose(f); + + + // Free D_inv + block_relax_free(D_inv); + + delete A; + delete A_sten; + + + +} // end of TEST(AnisoSpMVTest, TestsInUtil) // + diff --git a/test_data/RAPtorTests.ipynb b/test_data/RAPtorTests.ipynb index 7988cc29..bd40eba2 100644 --- a/test_data/RAPtorTests.ipynb +++ b/test_data/RAPtorTests.ipynb @@ -755,7 +755,7 @@ }, { "cell_type": "code", - "execution_count": 69, + "execution_count": 75, "metadata": {}, "outputs": [], "source": [