GRASS Programmer's Manual  6.4.2(2012)
test_solvers.c
Go to the documentation of this file.
00001 
00002 /*****************************************************************************
00003 *
00004 * MODULE:       Grass PDE Numerical Library
00005 * AUTHOR(S):    Soeren Gebbert, Berlin (GER) Dec 2006
00006 *               soerengebbert <at> gmx <dot> de
00007 *               
00008 * PURPOSE:      Unit tests for les solving
00009 *
00010 * COPYRIGHT:    (C) 2000 by the GRASS Development Team
00011 *
00012 *               This program is free software under the GNU General Public
00013 *               License (>=v2). Read the file COPYING that comes with GRASS
00014 *               for details.
00015 *
00016 *****************************************************************************/
00017 
00018 #include <stdio.h>
00019 #include <stdlib.h>
00020 #include <string.h>
00021 #include <grass/glocale.h>
00022 #include <grass/N_pde.h>
00023 #include "test_gpde_lib.h"
00024 
00025 /* prototypes */
00026 static int test_solvers(void);
00027 static N_les *create_normal_les(int rows);
00028 static N_les *create_sparse_les(int rows);
00029 
00030 /* ************************************************************************* */
00031 /* Performe the solver unit tests ****************************************** */
00032 /* ************************************************************************* */
00033 int unit_test_solvers(void)
00034 {
00035     int sum = 0;
00036 
00037     G_message(_("\n++ Running solver unit tests ++"));
00038 
00039     sum += test_solvers();
00040 
00041     if (sum > 0)
00042         G_warning(_("\n-- Solver unit tests failure --"));
00043     else
00044         G_message(_("\n-- Solver unit tests finished successfully --"));
00045 
00046     return sum;
00047 }
00048 
00049 /* *************************************************************** */
00050 /* create a normal matrix with values ** Hilbert matrix ********** */
00051 /* *************************************************************** */
00052 N_les *create_normal_les(int rows)
00053 {
00054     N_les *les;
00055     int i, j;
00056     int size = rows;
00057     double val;
00058 
00059     les = N_alloc_les(rows, N_NORMAL_LES);
00060     for (i = 0; i < size; i++) {
00061         val = 0.0;
00062         for (j = 0; j < size; j++) {
00063             les->A[i][j] = (double)(1.0 / (((double)i + 1.0) +
00064                                            ((double)j + 1.0) - 1.0));
00065             val += les->A[i][j];
00066         }
00067         les->b[i] = val;
00068     }
00069 
00070     return les;
00071 }
00072 
00073 /* *************************************************************** */
00074 /* create a sparse matrix with values ** Hilbert matrix ********** */
00075 /* *************************************************************** */
00076 N_les *create_sparse_les(int rows)
00077 {
00078     N_les *les;
00079     N_spvector *spvector;
00080     int i, j;
00081     double val;
00082 
00083     les = N_alloc_les(rows, N_SPARSE_LES);
00084 
00085     for (i = 0; i < rows; i++) {
00086         spvector = N_alloc_spvector(rows);
00087         val = 0;
00088 
00089         for (j = 0; j < rows; j++) {
00090             spvector->values[j] =
00091                 (double)(1.0 / (((double)i + 1.0) + ((double)j + 1.0) - 1.0));
00092             spvector->index[j] = j;
00093             val += spvector->values[j];
00094         }
00095 
00096         N_add_spvector_to_les(les, spvector, i);
00097         les->b[i] = val;
00098     }
00099 
00100 
00101     return les;
00102 }
00103 
00104 
00105 /* *************************************************************** */
00106 /* Test all implemented solvers for sparse and normal matrices *** */
00107 /* *************************************************************** */
00108 int test_solvers(void)
00109 {
00110     N_les *les;
00111     N_les *sples;
00112 
00113     G_message("\t * testing jacobi solver\n");
00114 
00115     les = create_normal_les(TEST_N_NUM_ROWS);
00116     sples = create_sparse_les(TEST_N_NUM_ROWS);
00117 
00118     N_solver_jacobi(les, 100, 1, 0.1e-4);
00119     /*N_print_les(les); */
00120     N_solver_jacobi(sples, 100, 1, 0.1e-4);
00121     /*N_print_les(sples); */
00122 
00123     N_free_les(les);
00124     N_free_les(sples);
00125 
00126 
00127     G_message("\t * testing SOR solver\n");
00128 
00129     les = create_normal_les(TEST_N_NUM_ROWS);
00130     sples = create_sparse_les(TEST_N_NUM_ROWS);
00131 
00132     N_solver_SOR(les, 100, 1, 0.1e-4);
00133     /*N_print_les(les); */
00134     N_solver_SOR(sples, 100, 1, 0.1e-4);
00135     /*N_print_les(sples); */
00136 
00137     N_free_les(les);
00138     N_free_les(sples);
00139 
00140     G_message("\t * testing cg solver\n");
00141 
00142     les = create_normal_les(TEST_N_NUM_ROWS);
00143     sples = create_sparse_les(TEST_N_NUM_ROWS);
00144 
00145     N_solver_cg(les, 100, 0.1e-8);
00146     /*N_print_les(les); */
00147     N_solver_cg(sples, 100, 0.1e-8);
00148     /*N_print_les(sples); */
00149 
00150     N_free_les(les);
00151     N_free_les(sples);
00152 
00153     G_message("\t * testing pcg solver with N_DIAGONAL_PRECONDITION\n");
00154 
00155     les = create_normal_les(TEST_N_NUM_ROWS);
00156     sples = create_sparse_les(TEST_N_NUM_ROWS);
00157 
00158     N_solver_pcg(les, 100, 0.1e-8, N_DIAGONAL_PRECONDITION);
00159     N_print_les(les);
00160     N_solver_pcg(sples, 100, 0.1e-8, N_DIAGONAL_PRECONDITION);
00161     N_print_les(sples);
00162 
00163     N_free_les(les);
00164     N_free_les(sples);
00165 
00166     G_message
00167         ("\t * testing pcg solver with N_ROWSCALE_EUKLIDNORM_PRECONDITION\n");
00168 
00169     les = create_normal_les(TEST_N_NUM_ROWS);
00170     sples = create_sparse_les(TEST_N_NUM_ROWS);
00171 
00172     N_solver_pcg(les, 100, 0.1e-8, N_ROWSCALE_EUKLIDNORM_PRECONDITION);
00173     N_print_les(les);
00174     N_solver_pcg(sples, 100, 0.1e-8, N_ROWSCALE_EUKLIDNORM_PRECONDITION);
00175     N_print_les(sples);
00176 
00177     N_free_les(les);
00178     N_free_les(sples);
00179 
00180     G_message
00181         ("\t * testing pcg solver with N_ROWSCALE_ABSSUMNORM_PRECONDITION\n");
00182 
00183     les = create_normal_les(TEST_N_NUM_ROWS);
00184     sples = create_sparse_les(TEST_N_NUM_ROWS);
00185 
00186     N_solver_pcg(les, 100, 0.1e-8, N_ROWSCALE_ABSSUMNORM_PRECONDITION);
00187     N_print_les(les);
00188     N_solver_pcg(sples, 100, 0.1e-8, N_ROWSCALE_ABSSUMNORM_PRECONDITION);
00189     N_print_les(sples);
00190 
00191     N_free_les(les);
00192     N_free_les(sples);
00193 
00194 
00195     G_message("\t * testing bicgstab solver\n");
00196 
00197     les = create_normal_les(TEST_N_NUM_ROWS);
00198     sples = create_sparse_les(TEST_N_NUM_ROWS);
00199 
00200     N_solver_bicgstab(les, 100, 0.1e-8);
00201     /*N_print_les(les); */
00202     N_solver_bicgstab(sples, 100, 0.1e-8);
00203     /*N_print_les(sples); */
00204 
00205     N_free_les(les);
00206     N_free_les(sples);
00207 
00208     G_message("\t * testing gauss elimination solver\n");
00209 
00210     les = create_normal_les(TEST_N_NUM_ROWS);
00211 
00212      /*GAUSS*/ N_solver_gauss(les);
00213     N_print_les(les);
00214 
00215     N_free_les(les);
00216 
00217     G_message("\t * testing lu decomposition solver\n");
00218 
00219     les = create_normal_les(TEST_N_NUM_ROWS);
00220 
00221      /*LU*/ N_solver_lu(les);
00222     N_print_les(les);
00223 
00224     N_free_les(les);
00225 
00226     les = create_normal_les(TEST_N_NUM_ROWS);
00227 
00228     /*cholesky */ N_solver_cholesky(les);
00229     N_print_les(les);
00230 
00231     N_free_les(les);
00232 
00233 
00234     return 0;
00235 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines