From 12a8f65ccb1aeb97aad60e1aab101ca9889a453d Mon Sep 17 00:00:00 2001 From: Amanda Bienz Date: Tue, 21 May 2024 21:05:22 -0600 Subject: [PATCH] Fixed GS tests --- raptor/util/linalg/relax.cpp | 26 ++---- raptor/util/tests/test_gs_aniso.cpp | 9 +- raptor/util/tests/test_gs_laplacian.cpp | 9 +- test_data/RAPtorTests.ipynb | 110 ++++++++++++++++++++++-- 4 files changed, 118 insertions(+), 36 deletions(-) diff --git a/raptor/util/linalg/relax.cpp b/raptor/util/linalg/relax.cpp index 1e2c450f..62159d8a 100644 --- a/raptor/util/linalg/relax.cpp +++ b/raptor/util/linalg/relax.cpp @@ -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; @@ -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; @@ -72,22 +66,19 @@ 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; } } } @@ -95,9 +86,6 @@ void sor(CSRMatrix* A, Vector& b, Vector& x, Vector& tmp, int num_sweeps, 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; diff --git a/raptor/util/tests/test_gs_aniso.cpp b/raptor/util/tests/test_gs_aniso.cpp index 8a85d035..93d28d5d 100644 --- a/raptor/util/tests/test_gs_aniso.cpp +++ b/raptor/util/tests/test_gs_aniso.cpp @@ -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 @@ -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); diff --git a/raptor/util/tests/test_gs_laplacian.cpp b/raptor/util/tests/test_gs_laplacian.cpp index 69619886..cd089578 100644 --- a/raptor/util/tests/test_gs_laplacian.cpp +++ b/raptor/util/tests/test_gs_laplacian.cpp @@ -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 @@ -89,6 +89,5 @@ TEST(AnisoJacobiTest, TestsInUtil) } fclose(f); - } // end of TEST(AnisoSpMVTest, TestsInUtil) // diff --git a/test_data/RAPtorTests.ipynb b/test_data/RAPtorTests.ipynb index 59a7de05..07e063c4 100644 --- a/test_data/RAPtorTests.ipynb +++ b/test_data/RAPtorTests.ipynb @@ -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" @@ -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", @@ -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", @@ -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",