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: linear equation system solvers 00009 * part of the gpde library 00010 * 00011 * COPYRIGHT: (C) 2000 by the GRASS Development Team 00012 * 00013 * This program is free software under the GNU General Public 00014 * License (>=v2). Read the file COPYING that comes with GRASS 00015 * for details. 00016 * 00017 *****************************************************************************/ 00018 00019 #include <math.h> 00020 #include <unistd.h> 00021 #include <stdio.h> 00022 #include <string.h> 00023 #include "grass/N_pde.h" 00024 #include "solvers_local_proto.h" 00025 00026 /*local protos */ 00027 static void scalar_product(double *a, double *b, double *scalar, int rows); 00028 static void sub_vectors(double *source_a, double *source_b, double *result, 00029 int row); 00030 static void sub_vectors_scalar(double *source_a, double *source_b, 00031 double *result, double scalar_b, int rows); 00032 static void add_vectors(double *source_a, double *source_b, double *result, 00033 int rows); 00034 static void add_vectors_scalar(double *source_a, double *source_b, 00035 double *result, double scalar_b, int rows); 00036 static void add_vectors_scalar2(double *source_a, double *source_b, 00037 double *result, double scalar_a, 00038 double scalar_b, int rows); 00039 static void scalar_vector_product(double *a, double *result, double scalar, 00040 int rows); 00041 static void sync_vectors(double *source, double *target, int rows); 00042 00043 00044 /* ******************************************************* * 00045 * *** preconditioned conjugate gradients **************** * 00046 * ******************************************************* */ 00064 int N_solver_pcg(N_les * L, int maxit, double err, int prec) 00065 { 00066 double *r, *z; 00067 double *p; 00068 double *v; 00069 double *x, *b; 00070 double s = 0.0; 00071 double a0 = 0, a1 = 0, mygamma, tmp = 0; 00072 int m, rows, i; 00073 int finished = 2; 00074 int error_break; 00075 N_les *M; 00076 00077 if (L->quad != 1) { 00078 G_warning(_("The linear equation system is not quadratic")); 00079 return -1; 00080 } 00081 00082 /* check for symmetry */ 00083 if (check_symmetry(L) != 1) { 00084 G_warning(_("Matrix is not symmetric!")); 00085 } 00086 00087 x = L->x; 00088 b = L->b; 00089 rows = L->rows; 00090 00091 r = vectmem(rows); 00092 p = vectmem(rows); 00093 v = vectmem(rows); 00094 z = vectmem(rows); 00095 00096 error_break = 0; 00097 00098 /*compute the preconditioning matrix */ 00099 M = N_create_diag_precond_matrix(L, prec); 00100 00101 /* 00102 * residual calculation 00103 */ 00104 #pragma omp parallel 00105 { 00106 /* matrix vector multiplication */ 00107 if (L->type == N_SPARSE_LES) 00108 N_sparse_matrix_vector_product(L, x, v); 00109 else 00110 N_matrix_vector_product(L, x, v); 00111 00112 sub_vectors(b, v, r, rows); 00113 /*performe the preconditioning */ 00114 N_sparse_matrix_vector_product(M, r, p); 00115 00116 /* scalar product */ 00117 #pragma omp for schedule (static) private(i) reduction(+:s) 00118 for (i = 0; i < rows; i++) { 00119 s += p[i] * r[i]; 00120 } 00121 } 00122 00123 a0 = s; 00124 s = 0.0; 00125 00126 /* ******************* */ 00127 /* start the iteration */ 00128 /* ******************* */ 00129 for (m = 0; m < maxit; m++) { 00130 #pragma omp parallel default(shared) 00131 { 00132 /* matrix vector multiplication */ 00133 if (L->type == N_SPARSE_LES) 00134 N_sparse_matrix_vector_product(L, p, v); 00135 else 00136 N_matrix_vector_product(L, p, v); 00137 00138 00139 /* scalar product */ 00140 #pragma omp for schedule (static) private(i) reduction(+:s) 00141 for (i = 0; i < rows; i++) { 00142 s += v[i] * p[i]; 00143 } 00144 00145 /* barrier */ 00146 #pragma omp single 00147 { 00148 tmp = s; 00149 mygamma = a0 / tmp; 00150 s = 0.0; 00151 } 00152 00153 add_vectors_scalar(x, p, x, mygamma, rows); 00154 00155 if (m % 50 == 1) { 00156 /* matrix vector multiplication */ 00157 if (L->type == N_SPARSE_LES) 00158 N_sparse_matrix_vector_product(L, x, v); 00159 else 00160 N_matrix_vector_product(L, x, v); 00161 00162 sub_vectors(b, v, r, rows); 00163 00164 } 00165 else { 00166 sub_vectors_scalar(r, v, r, mygamma, rows); 00167 } 00168 00169 /*performe the preconditioning */ 00170 N_sparse_matrix_vector_product(M, r, z); 00171 00172 /* scalar product */ 00173 #pragma omp for schedule (static) private(i) reduction(+:s) 00174 for (i = 0; i < rows; i++) { 00175 s += z[i] * r[i]; 00176 } 00177 00178 /* barrier */ 00179 #pragma omp single 00180 { 00181 a1 = s; 00182 tmp = a1 / a0; 00183 a0 = a1; 00184 s = 0.0; 00185 00186 if (a1 < 0 || a1 == 0 || a1 > 0) { 00187 ; 00188 } 00189 else { 00190 G_warning(_("Unable to solve the linear equation system")); 00191 error_break = 1; 00192 } 00193 } 00194 add_vectors_scalar(z, p, p, tmp, rows); 00195 } 00196 00197 if (L->type == N_SPARSE_LES) 00198 G_message(_("Sparse PCG -- iteration %i error %g\n"), m, a0); 00199 else 00200 G_message(_("PCG -- iteration %i error %g\n"), m, a0); 00201 00202 if (error_break == 1) { 00203 finished = -1; 00204 break; 00205 } 00206 00207 00208 if (a0 < err) { 00209 finished = 1; 00210 break; 00211 } 00212 } 00213 00214 G_free(r); 00215 G_free(p); 00216 G_free(v); 00217 G_free(z); 00218 00219 return finished; 00220 } 00221 00222 00223 /* ******************************************************* * 00224 * ****************** conjugate gradients **************** * 00225 * ******************************************************* */ 00242 int N_solver_cg(N_les * L, int maxit, double err) 00243 { 00244 double *r; 00245 double *p; 00246 double *v; 00247 double *x, *b; 00248 double s = 0.0; 00249 double a0 = 0, a1 = 0, mygamma, tmp = 0; 00250 int m, rows, i; 00251 int finished = 2; 00252 int error_break; 00253 00254 if (L->quad != 1) { 00255 G_warning(_("The linear equation system is not quadratic")); 00256 return -1; 00257 } 00258 00259 /* check for symmetry */ 00260 if (check_symmetry(L) != 1) { 00261 G_warning(_("Matrix is not symmetric!")); 00262 } 00263 00264 x = L->x; 00265 b = L->b; 00266 rows = L->rows; 00267 00268 r = vectmem(rows); 00269 p = vectmem(rows); 00270 v = vectmem(rows); 00271 00272 error_break = 0; 00273 /* 00274 * residual calculation 00275 */ 00276 #pragma omp parallel 00277 { 00278 /* matrix vector multiplication */ 00279 if (L->type == N_SPARSE_LES) 00280 N_sparse_matrix_vector_product(L, x, v); 00281 else 00282 N_matrix_vector_product(L, x, v); 00283 00284 sub_vectors(b, v, r, rows); 00285 sync_vectors(r, p, rows); 00286 00287 /* scalar product */ 00288 #pragma omp for schedule (static) private(i) reduction(+:s) 00289 for (i = 0; i < rows; i++) { 00290 s += r[i] * r[i]; 00291 } 00292 } 00293 00294 a0 = s; 00295 s = 0.0; 00296 00297 /* ******************* */ 00298 /* start the iteration */ 00299 /* ******************* */ 00300 for (m = 0; m < maxit; m++) { 00301 #pragma omp parallel default(shared) 00302 { 00303 /* matrix vector multiplication */ 00304 if (L->type == N_SPARSE_LES) 00305 N_sparse_matrix_vector_product(L, p, v); 00306 else 00307 N_matrix_vector_product(L, p, v); 00308 00309 00310 /* scalar product */ 00311 #pragma omp for schedule (static) private(i) reduction(+:s) 00312 for (i = 0; i < rows; i++) { 00313 s += v[i] * p[i]; 00314 } 00315 00316 /* barrier */ 00317 #pragma omp single 00318 { 00319 tmp = s; 00320 mygamma = a0 / tmp; 00321 s = 0.0; 00322 } 00323 00324 add_vectors_scalar(x, p, x, mygamma, rows); 00325 00326 if (m % 50 == 1) { 00327 /* matrix vector multiplication */ 00328 if (L->type == N_SPARSE_LES) 00329 N_sparse_matrix_vector_product(L, x, v); 00330 else 00331 N_matrix_vector_product(L, x, v); 00332 00333 sub_vectors(b, v, r, rows); 00334 } 00335 else { 00336 sub_vectors_scalar(r, v, r, mygamma, rows); 00337 } 00338 00339 /* scalar product */ 00340 #pragma omp for schedule (static) private(i) reduction(+:s) 00341 for (i = 0; i < rows; i++) { 00342 s += r[i] * r[i]; 00343 } 00344 00345 /* barrier */ 00346 #pragma omp single 00347 { 00348 a1 = s; 00349 tmp = a1 / a0; 00350 a0 = a1; 00351 s = 0.0; 00352 00353 if (a1 < 0 || a1 == 0 || a1 > 0) { 00354 ; 00355 } 00356 else { 00357 G_warning(_("Unable to solve the linear equation system")); 00358 error_break = 1; 00359 } 00360 } 00361 add_vectors_scalar(r, p, p, tmp, rows); 00362 } 00363 00364 if (L->type == N_SPARSE_LES) 00365 G_message(_("Sparse CG -- iteration %i error %g\n"), m, a0); 00366 else 00367 G_message(_("CG -- iteration %i error %g\n"), m, a0); 00368 00369 if (error_break == 1) { 00370 finished = -1; 00371 break; 00372 } 00373 00374 if (a0 < err) { 00375 finished = 1; 00376 break; 00377 } 00378 } 00379 00380 G_free(r); 00381 G_free(p); 00382 G_free(v); 00383 00384 return finished; 00385 } 00386 00387 /* ******************************************************* * 00388 * ************ biconjugate gradients ******************** * 00389 * ******************************************************* */ 00407 int N_solver_bicgstab(N_les * L, int maxit, double err) 00408 { 00409 double *r; 00410 double *r0; 00411 double *p; 00412 double *v; 00413 double *s; 00414 double *t; 00415 double *x, *b; 00416 double s1 = 0.0, s2 = 0.0, s3 = 0.0; 00417 double alpha = 0, beta = 0, omega, rr0 = 0, error; 00418 int m, rows, i; 00419 int finished = 2; 00420 int error_break; 00421 00422 if (L->quad != 1) { 00423 G_warning(_("The linear equation system is not quadratic")); 00424 return -1; 00425 } 00426 00427 00428 x = L->x; 00429 b = L->b; 00430 rows = L->rows; 00431 r = vectmem(rows); 00432 r0 = vectmem(rows); 00433 p = vectmem(rows); 00434 v = vectmem(rows); 00435 s = vectmem(rows); 00436 t = vectmem(rows); 00437 00438 error_break = 0; 00439 00440 #pragma omp parallel 00441 { 00442 if (L->type == N_SPARSE_LES) 00443 N_sparse_matrix_vector_product(L, x, v); 00444 else 00445 N_matrix_vector_product(L, x, v); 00446 sub_vectors(b, v, r, rows); 00447 sync_vectors(r, r0, rows); 00448 sync_vectors(r, p, rows); 00449 } 00450 00451 s1 = s2 = s3 = 0.0; 00452 00453 /* ******************* */ 00454 /* start the iteration */ 00455 /* ******************* */ 00456 for (m = 0; m < maxit; m++) { 00457 00458 #pragma omp parallel default(shared) 00459 { 00460 if (L->type == N_SPARSE_LES) 00461 N_sparse_matrix_vector_product(L, p, v); 00462 else 00463 N_matrix_vector_product(L, p, v); 00464 00465 /* scalar product */ 00466 #pragma omp for schedule (static) private(i) reduction(+:s1, s2, s3) 00467 for (i = 0; i < rows; i++) { 00468 s1 += r[i] * r[i]; 00469 s2 += r[i] * r0[i]; 00470 s3 += v[i] * r0[i]; 00471 } 00472 00473 #pragma omp single 00474 { 00475 error = s1; 00476 00477 if (error < 0 || error == 0 || error > 0) { 00478 ; 00479 } 00480 else { 00481 G_warning(_("Unable to solve the linear equation system")); 00482 error_break = 1; 00483 } 00484 00485 rr0 = s2; 00486 alpha = rr0 / s3; 00487 s1 = s2 = s3 = 0.0; 00488 } 00489 00490 sub_vectors_scalar(r, v, s, alpha, rows); 00491 if (L->type == N_SPARSE_LES) 00492 N_sparse_matrix_vector_product(L, s, t); 00493 else 00494 N_matrix_vector_product(L, s, t); 00495 00496 /* scalar product */ 00497 #pragma omp for schedule (static) private(i) reduction(+:s1, s2) 00498 for (i = 0; i < rows; i++) { 00499 s1 += t[i] * s[i]; 00500 s2 += t[i] * t[i]; 00501 } 00502 00503 #pragma omp single 00504 { 00505 omega = s1 / s2; 00506 s1 = s2 = 0.0; 00507 } 00508 00509 add_vectors_scalar2(p, s, r, alpha, omega, rows); 00510 add_vectors(x, r, x, rows); 00511 sub_vectors_scalar(s, t, r, omega, rows); 00512 00513 #pragma omp for schedule (static) private(i) reduction(+:s1) 00514 for (i = 0; i < rows; i++) { 00515 s1 += r[i] * r0[i]; 00516 } 00517 00518 #pragma omp single 00519 { 00520 beta = alpha / omega * s1 / rr0; 00521 s1 = s2 = s3 = 0.0; 00522 } 00523 00524 sub_vectors_scalar(p, v, p, omega, rows); 00525 add_vectors_scalar(r, p, p, beta, rows); 00526 } 00527 00528 00529 if (L->type == N_SPARSE_LES) 00530 G_message(_("Sparse BiCGStab -- iteration %i error %g\n"), m, 00531 error); 00532 else 00533 G_message(_("BiCGStab -- iteration %i error %g\n"), m, error); 00534 00535 if (error_break == 1) { 00536 finished = -1; 00537 break; 00538 } 00539 00540 if (error < err) { 00541 finished = 1; 00542 break; 00543 } 00544 } 00545 00546 G_free(r); 00547 G_free(r0); 00548 G_free(p); 00549 G_free(v); 00550 G_free(s); 00551 G_free(t); 00552 00553 return finished; 00554 } 00555 00568 void scalar_product(double *a, double *b, double *scalar, int rows) 00569 { 00570 int i; 00571 double s = 0.0; 00572 00573 #pragma omp parallel for schedule (static) reduction(+:s) 00574 for (i = 0; i < rows; i++) { 00575 s += a[i] * b[i]; 00576 } 00577 00578 *scalar = s; 00579 return; 00580 } 00581 00594 void N_matrix_vector_product(N_les * L, double *x, double *result) 00595 { 00596 int i, j; 00597 double tmp; 00598 00599 #pragma omp for schedule (static) private(i, j, tmp) 00600 for (i = 0; i < L->rows; i++) { 00601 tmp = 0; 00602 for (j = 0; j < L->cols; j++) { 00603 tmp += L->A[i][j] * x[j]; 00604 } 00605 result[i] = tmp; 00606 } 00607 return; 00608 } 00609 00622 void N_sparse_matrix_vector_product(N_les * L, double *x, double *result) 00623 { 00624 int i, j; 00625 double tmp; 00626 00627 #pragma omp for schedule (static) private(i, j, tmp) 00628 for (i = 0; i < L->rows; i++) { 00629 tmp = 0; 00630 for (j = 0; j < L->Asp[i]->cols; j++) { 00631 tmp += L->Asp[i]->values[j] * x[L->Asp[i]->index[j]]; 00632 } 00633 result[i] = tmp; 00634 } 00635 return; 00636 } 00637 00652 void 00653 add_vectors_scalar2(double *a, double *b, double *result, double scalar_a, 00654 double scalar_b, int rows) 00655 { 00656 int i; 00657 00658 #pragma omp for schedule (static) 00659 for (i = 0; i < rows; i++) { 00660 result[i] = scalar_a * a[i] + scalar_b * b[i]; 00661 } 00662 00663 return; 00664 } 00665 00679 void 00680 add_vectors_scalar(double *a, double *b, double *result, double scalar_b, 00681 int rows) 00682 { 00683 int i; 00684 00685 #pragma omp for schedule (static) 00686 for (i = 0; i < rows; i++) { 00687 result[i] = a[i] + scalar_b * b[i]; 00688 } 00689 00690 return; 00691 } 00692 00706 void 00707 sub_vectors_scalar(double *a, double *b, double *result, double scalar_b, 00708 int rows) 00709 { 00710 int i; 00711 00712 #pragma omp for schedule (static) 00713 for (i = 0; i < rows; i++) { 00714 result[i] = a[i] - scalar_b * b[i]; 00715 } 00716 00717 return; 00718 } 00719 00732 void add_vectors(double *a, double *b, double *result, int rows) 00733 { 00734 int i; 00735 00736 #pragma omp for schedule (static) 00737 for (i = 0; i < rows; i++) { 00738 result[i] = a[i] + b[i]; 00739 } 00740 00741 return; 00742 } 00743 00756 void sub_vectors(double *a, double *b, double *result, int rows) 00757 { 00758 int i; 00759 00760 #pragma omp for schedule (static) 00761 for (i = 0; i < rows; i++) { 00762 result[i] = a[i] - b[i]; 00763 } 00764 00765 return; 00766 } 00767 00780 void scalar_vector_product(double *a, double *result, double scalar, int rows) 00781 { 00782 int i; 00783 00784 #pragma omp for schedule (static) 00785 for (i = 0; i < rows; i++) { 00786 result[i] = scalar * a[i]; 00787 } 00788 00789 return; 00790 } 00791 00800 void sync_vectors(double *source, double *target, int rows) 00801 { 00802 int i; 00803 00804 #pragma omp for schedule (static) 00805 for (i = 0; i < rows; i++) { 00806 target[i] = source[i]; 00807 } 00808 00809 return; 00810 } 00811 00812 /* ******************************************************* * 00813 * ***** Check if matrix is symmetric ******************** * 00814 * ******************************************************* */ 00823 int check_symmetry(N_les * L) 00824 { 00825 int i, j, k; 00826 double value1 = 0; 00827 double value2 = 0; 00828 int index; 00829 int symm = 0; 00830 00831 if (L->quad != 1) { 00832 G_warning(_("The linear equation system is not quadratic")); 00833 return 0; 00834 } 00835 00836 G_debug(2, "check_symmetry: Check if matrix is symmetric"); 00837 00838 if (L->type == N_SPARSE_LES) { 00839 #pragma omp parallel for schedule (static) private(i, j, k, value1, value2, index) reduction(+:symm) shared(L) 00840 for (j = 0; j < L->rows; j++) { 00841 for (i = 1; i < L->Asp[j]->cols; i++) { 00842 value1 = 0; 00843 value2 = 0; 00844 index = 0; 00845 00846 index = L->Asp[j]->index[i]; 00847 value1 = L->Asp[j]->values[i]; 00848 00849 for (k = 1; k < L->Asp[index]->cols; k++) { 00850 if (L->Asp[index]->index[k] == j) { 00851 value2 = L->Asp[index]->values[k]; 00852 if ((value1 != value2)) { 00853 if ((fabs((fabs(value1) - fabs(value2))) < 00854 SYMM_TOLERANCE)) { 00855 G_debug(5, 00856 "check_symmetry: sparse matrix is unsymmetric, but within tolerance"); 00857 } 00858 else { 00859 G_warning 00860 ("Matrix unsymmetric: Position [%i][%i] : [%i][%i] \nError: %12.18lf != %12.18lf \ndifference = %12.18lf\nStop symmetry calculation.\n", 00861 j, index, index, L->Asp[index]->index[k], 00862 value1, value2, 00863 fabs(fabs(value1) - fabs(value2))); 00864 symm++; 00865 } 00866 } 00867 } 00868 } 00869 } 00870 } 00871 } 00872 else { 00873 #pragma omp parallel for schedule (static) private(i, j, k, value1, value2, index) reduction(+:symm) shared(L) 00874 for (i = 0; i < L->rows; i++) { 00875 for (j = i + 1; j < L->rows; j++) { 00876 if (L->A[i][j] != L->A[j][i]) { 00877 if ((fabs(fabs(L->A[i][j]) - fabs(L->A[j][i])) < 00878 SYMM_TOLERANCE)) { 00879 G_debug(5, 00880 "check_symmetry: matrix is unsymmetric, but within tolerance"); 00881 } 00882 else { 00883 G_warning 00884 ("Matrix unsymmetric: Position [%i][%i] : [%i][%i] \nError: %12.18lf != %12.18lf\ndifference = %12.18lf\nStop symmetry calculation.\n", 00885 i, j, j, i, L->A[i][j], L->A[j][i], 00886 fabs(fabs(L->A[i][j]) - fabs(L->A[j][i]))); 00887 symm++; 00888 } 00889 } 00890 } 00891 } 00892 } 00893 00894 if (symm > 0) 00895 return 0; 00896 00897 return 1; 00898 } 00899 00900 00909 N_les *N_create_diag_precond_matrix(N_les * L, int prec) 00910 { 00911 N_les *L_new; 00912 int rows = L->rows; 00913 int cols = L->cols; 00914 int i, j; 00915 double sum; 00916 00917 L_new = N_alloc_les_A(rows, N_SPARSE_LES); 00918 00919 if (L->type == N_NORMAL_LES) { 00920 #pragma omp parallel for schedule (static) private(i, sum) shared(L_new, L, rows, prec) 00921 for (i = 0; i < rows; i++) { 00922 N_spvector *spvect = N_alloc_spvector(1); 00923 00924 switch (prec) { 00925 case N_ROWSCALE_EUKLIDNORM_PRECONDITION: 00926 sum = 0; 00927 for (j = 0; j < cols; j++) 00928 sum += L->A[i][j] * L->A[i][j]; 00929 spvect->values[0] = 1.0 / sqrt(sum); 00930 break; 00931 case N_ROWSCALE_ABSSUMNORM_PRECONDITION: 00932 sum = 0; 00933 for (j = 0; j < cols; j++) 00934 sum += fabs(L->A[i][j]); 00935 spvect->values[0] = 1.0 / (sum); 00936 break; 00937 case N_DIAGONAL_PRECONDITION: 00938 spvect->values[0] = 1.0 / L->A[i][i]; 00939 break; 00940 default: 00941 spvect->values[0] = 1.0 / L->A[i][i]; 00942 } 00943 00944 00945 spvect->index[0] = i; 00946 spvect->cols = 1;; 00947 N_add_spvector_to_les(L_new, spvect, i); 00948 00949 } 00950 } 00951 else { 00952 #pragma omp parallel for schedule (static) private(i, sum) shared(L_new, L, rows, prec) 00953 for (i = 0; i < rows; i++) { 00954 N_spvector *spvect = N_alloc_spvector(1); 00955 00956 switch (prec) { 00957 case N_ROWSCALE_EUKLIDNORM_PRECONDITION: 00958 sum = 0; 00959 for (j = 0; j < L->Asp[i]->cols; j++) 00960 sum += L->Asp[i]->values[j] * L->Asp[i]->values[j]; 00961 spvect->values[0] = 1.0 / sqrt(sum); 00962 break; 00963 case N_ROWSCALE_ABSSUMNORM_PRECONDITION: 00964 sum = 0; 00965 for (j = 0; j < L->Asp[i]->cols; j++) 00966 sum += fabs(L->Asp[i]->values[j]); 00967 spvect->values[0] = 1.0 / (sum); 00968 break; 00969 case N_DIAGONAL_PRECONDITION: 00970 spvect->values[0] = 1.0 / L->Asp[i]->values[0]; 00971 break; 00972 default: 00973 spvect->values[0] = 1.0 / L->Asp[i]->values[0]; 00974 } 00975 00976 spvect->index[0] = i; 00977 spvect->cols = 1;; 00978 N_add_spvector_to_les(L_new, spvect, i); 00979 } 00980 } 00981 return L_new; 00982 }