GRASS Programmer's Manual
6.4.2(2012)
|
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 }