GRASS Programmer's Manual  6.4.2(2012)
N_solvers_krylov.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:      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 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines