Skip to content

Commit

Permalink
Fixed GS tests
Browse files Browse the repository at this point in the history
  • Loading branch information
Amanda Bienz authored and Amanda Bienz committed May 22, 2024
1 parent 6ed68ed commit 12a8f65
Show file tree
Hide file tree
Showing 4 changed files with 118 additions and 36 deletions.
26 changes: 7 additions & 19 deletions raptor/util/linalg/relax.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,9 +19,6 @@ namespace raptor {
void jacobi(CSRMatrix* A, Vector& b, Vector& x, Vector& tmp, int num_sweeps,
double omega)
{
A->sort();
A->move_diag();

int row_start, row_end;
double diag, row_sum;

Expand Down Expand Up @@ -57,9 +54,6 @@ void jacobi(CSRMatrix* A, Vector& b, Vector& x, Vector& tmp, int num_sweeps,
void sor(CSRMatrix* A, Vector& b, Vector& x, Vector& tmp, int num_sweeps,
double omega)
{
A->sort();
A->move_diag();

int row_start, row_end;
double diag;
double rsum;
Expand All @@ -72,32 +66,26 @@ void sor(CSRMatrix* A, Vector& b, Vector& x, Vector& tmp, int num_sweeps,
diag = 0;
row_start = A->idx1[i];
row_end = A->idx1[i+1];
if (row_start == row_end) continue;

if (A->idx2[row_start] == i)
{
diag = A->vals[row_start];
row_start;
}
else continue;

for (int j = row_start; j < row_end; j++)
{
rsum += A->vals[j] * x[A->idx2[j]];
int col = A->idx2[j];
if (i == col)
diag = A->vals[j];
else
rsum += A->vals[j] * x[col];
}

if (diag)
if (fabs(diag) > zero_tol)
x[i] = omega*(b[i] - rsum)/diag + (1 - omega) * x[i];
// x[i] = (b[i] - rsum) / diag;
}
}
}

void ssor(CSRMatrix* A, Vector& b, Vector& x, Vector& tmp, int num_sweeps,
double omega)
{
A->sort();
A->move_diag();

int row_start, row_end;
double diag_inv;
double orig_x = 0;
Expand Down
9 changes: 4 additions & 5 deletions raptor/util/tests/test_gs_aniso.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,10 +29,10 @@ TEST(AnisoJacobiTest, TestsInUtil)
Vector b(A_sten->n_rows);
Vector tmp(A_sten->n_rows);

const char* x_ones_1 = "../../../../test_data/aniso_sor_ones_1.txt";
const char* x_ones_2 = "../../../../test_data/aniso_sor_ones_2.txt";
const char* x_inc_1 = "../../../../test_data/aniso_sor_inc_1.txt";
const char* x_inc_2 = "../../../../test_data/aniso_sor_inc_2.txt";
const char* x_ones_1 = "../../../../test_data/aniso_gs_ones_1.txt";
const char* x_ones_2 = "../../../../test_data/aniso_gs_ones_2.txt";
const char* x_inc_1 = "../../../../test_data/aniso_gs_inc_1.txt";
const char* x_inc_2 = "../../../../test_data/aniso_gs_inc_2.txt";

/*********************************************
* Test sor when b is constant value 1.0
Expand All @@ -44,7 +44,6 @@ TEST(AnisoJacobiTest, TestsInUtil)
f = fopen(x_ones_1, "r");
for (int i = 0; i < A_sten->n_rows; i++)
{
printf("x[%d] = %f\n", i, x[i]);
n_items_read = fscanf(f, "%lg\n", &x_val);
ASSERT_EQ(n_items_read, 1);
ASSERT_NEAR(x[i],x_val,1e-06);
Expand Down
9 changes: 4 additions & 5 deletions raptor/util/tests/test_gs_laplacian.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,10 +27,10 @@ TEST(AnisoJacobiTest, TestsInUtil)
Vector b(A_sten->n_rows);
Vector tmp(A_sten->n_rows);

const char* x_ones_1 = "../../../../test_data/laplace_sor_ones_1.txt";
const char* x_ones_2 = "../../../../test_data/laplace_sor_ones_2.txt";
const char* x_inc_1 = "../../../../test_data/laplace_sor_inc_1.txt";
const char* x_inc_2 = "../../../../test_data/laplace_sor_inc_2.txt";
const char* x_ones_1 = "../../../../test_data/laplace_gs_ones_1.txt";
const char* x_ones_2 = "../../../../test_data/laplace_gs_ones_2.txt";
const char* x_inc_1 = "../../../../test_data/laplace_gs_inc_1.txt";
const char* x_inc_2 = "../../../../test_data/laplace_gs_inc_2.txt";

/*********************************************
* Test sor when b is constant value 1.0
Expand Down Expand Up @@ -89,6 +89,5 @@ TEST(AnisoJacobiTest, TestsInUtil)
}
fclose(f);


} // end of TEST(AnisoSpMVTest, TestsInUtil) //

110 changes: 103 additions & 7 deletions test_data/RAPtorTests.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,14 +2,14 @@
"cells": [
{
"cell_type": "code",
"execution_count": 2,
"execution_count": 49,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from pyamg.gallery import diffusion_stencil_2d, stencil_grid\n",
"from pyamg.classical.split import RS, PMIS, MIS\n",
"from pyamg.relaxation.relaxation import jacobi, block_jacobi, sor\n",
"from pyamg.relaxation.relaxation import jacobi, block_jacobi, sor, gauss_seidel\n",
"from scipy.io import mmwrite\n",
"from scipy.sparse import csr_matrix\n",
"import PetscBinaryIO"
Expand Down Expand Up @@ -755,9 +755,98 @@
},
{
"cell_type": "code",
"execution_count": 45,
"execution_count": 68,
"metadata": {},
"outputs": [],
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Val[0][0] = 2.600000e+01\n",
"\n",
"Val[0][1] = -1.000000e+00\n",
"\n",
"Val[0][2] = 0.000000e+00\n",
"\n",
"Val[0][3] = 0.000000e+00\n",
"\n",
"Val[0][4] = 0.000000e+00\n",
"\n",
"Val[0][10] = -1.000000e+00\n",
"\n",
"Val[0][11] = -1.000000e+00\n",
"\n",
"Val[0][12] = 0.000000e+00\n",
"\n",
"Val[0][13] = 0.000000e+00\n",
"\n",
"Val[0][14] = 0.000000e+00\n",
"\n",
"Val[0][100] = -1.000000e+00\n",
"\n",
"Val[0][101] = -1.000000e+00\n",
"\n",
"Val[0][102] = 0.000000e+00\n",
"\n",
"Val[0][103] = 0.000000e+00\n",
"\n",
"Val[0][104] = 0.000000e+00\n",
"\n",
"Val[0][110] = -1.000000e+00\n",
"\n",
"Val[0][111] = -1.000000e+00\n",
"\n",
"Val[0][112] = 0.000000e+00\n",
"\n",
"Val[0][113] = 0.000000e+00\n",
"\n",
"Val[0][114] = 0.000000e+00\n",
"\n",
"Val[0][5] = 0.000000e+00\n",
"\n",
"Val[0][6] = 0.000000e+00\n",
"\n",
"Val[0][7] = 0.000000e+00\n",
"\n",
"Val[0][8] = 0.000000e+00\n",
"\n",
"Val[0][9] = 0.000000e+00\n",
"\n",
"Val[0][15] = 0.000000e+00\n",
"\n",
"Val[0][16] = 0.000000e+00\n",
"\n",
"Val[0][17] = 0.000000e+00\n",
"\n",
"Val[0][18] = 0.000000e+00\n",
"\n",
"Val[0][19] = 0.000000e+00\n",
"\n",
"Val[0][105] = 0.000000e+00\n",
"\n",
"Val[0][106] = 0.000000e+00\n",
"\n",
"Val[0][107] = 0.000000e+00\n",
"\n",
"Val[0][108] = 0.000000e+00\n",
"\n",
"Val[0][109] = 0.000000e+00\n",
"\n",
"Val[0][115] = 0.000000e+00\n",
"\n",
"Val[0][116] = 0.000000e+00\n",
"\n",
"Val[0][117] = 0.000000e+00\n",
"\n",
"Val[0][118] = 0.000000e+00\n",
"\n",
"Val[0][119] = 0.000000e+00\n",
"\n",
"[1.]\n",
"[0.03846154]\n"
]
}
],
"source": [
"##############################\n",
"### Relaxation Tests\n",
Expand Down Expand Up @@ -795,12 +884,13 @@
"np.savetxt(\"laplace_jacobi_inc_2.txt\", x, fmt=\"%f\")\n",
"\n",
"## Forward GS\n",
"from scipy import sparse\n",
"Aniso = Aniso.tocsr()\n",
"b = np.ones((Aniso.shape[0],1))\n",
"x = np.zeros((Aniso.shape[0],1))\n",
"sor(Aniso, x, b, 1.0, iterations=1, sweep='forward')\n",
"gauss_seidel(Aniso, x, b)\n",
"np.savetxt(\"aniso_gs_ones_1.txt\", x, fmt=\"%f\")\n",
"sor(Aniso, x, b, 1.0, iterations=1, sweep='forward')\n",
"gauss_seidel(Aniso, x, b)\n",
"np.savetxt(\"aniso_gs_ones_2.txt\", x, fmt=\"%f\")\n",
"\n",
"b = np.arange(Aniso.shape[0], dtype='float')\n",
Expand All @@ -813,7 +903,13 @@
"Laplacian = Laplacian.tocsr()\n",
"b = np.ones((Laplacian.shape[0],1))\n",
"x = np.zeros((Laplacian.shape[0],1))\n",
"sor(Laplacian, x, b, 1.0, iterations=1, sweep='forward')\n",
"for i in range(Laplacian.shape[0]):\n",
" for j in range(Laplacian.indptr[i], Laplacian.indptr[i+1]):\n",
" print(\"Val[%d][%d] = %e\\n\" %(i, Laplacian.indices[j], Laplacian.data[j]))\n",
" break\n",
"print(b[0])\n",
"gauss_seidel(Laplacian, x, b)\n",
"print(x[0])\n",
"np.savetxt(\"laplace_gs_ones_1.txt\", x, fmt=\"%f\")\n",
"sor(Laplacian, x, b, 1.0, iterations=1, sweep='forward')\n",
"np.savetxt(\"laplace_gs_ones_2.txt\", x, fmt=\"%f\")\n",
Expand Down

0 comments on commit 12a8f65

Please sign in to comment.