Skip to content

Commit

Permalink
Added block jacobi tests
Browse files Browse the repository at this point in the history
  • Loading branch information
Amanda Bienz authored and Amanda Bienz committed May 24, 2024
1 parent 1c3457f commit 6e1cfba
Show file tree
Hide file tree
Showing 6 changed files with 164 additions and 30 deletions.
43 changes: 24 additions & 19 deletions raptor/util/linalg/relax.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand All @@ -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;
}
}
Expand All @@ -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;
Expand All @@ -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++)
Expand All @@ -205,33 +209,33 @@ 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++)
{
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]);

}
}
Expand Down Expand Up @@ -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];


}
}

Expand Down
4 changes: 2 additions & 2 deletions raptor/util/linalg/relax.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);


}
Expand Down
8 changes: 8 additions & 0 deletions raptor/util/tests/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
29 changes: 21 additions & 8 deletions raptor/util/tests/test_bsr_jacobi_aniso.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,18 +29,24 @@ 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
*********************************************/
// 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++)
{
Expand All @@ -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++)
{
Expand All @@ -60,6 +66,7 @@ TEST(AnisoJacobiTest, TestsInUtil)
ASSERT_NEAR(x[i],x_val,1e-06);
}
fclose(f);


/*********************************************
* Test jacobi when b[i] = i
Expand All @@ -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++)
{
Expand All @@ -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++)
{
Expand All @@ -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) //

108 changes: 108 additions & 0 deletions raptor/util/tests/test_bsr_jacobi_laplacian.cpp
View file
Original file line number Diff line number Diff line change
@@ -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) //

2 changes: 1 addition & 1 deletion test_data/RAPtorTests.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -755,7 +755,7 @@
},
{
"cell_type": "code",
"execution_count": 69,
"execution_count": 75,
"metadata": {},
"outputs": [],
"source": [
Expand Down

0 comments on commit 6e1cfba

Please sign in to comment.