Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
smsrres.c
Go to the documentation of this file.
1/*@HEADER
2// ***********************************************************************
3//
4// Ifpack: Object-Oriented Algebraic Preconditioner Package
5// Copyright (2002) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40//@HEADER
41*/
42
43#include <stdlib.h>
44#include <stdio.h>
45#include <math.h>
46#define max(x,y) (( x > y ) ? x : y) /* max function */
47double smsrres (int m, int n,
48 double *val, int *indx,
49 double *xlocal, double *x, double *b)
50{
51 int i, j, jbgn, jend, ione = 1;
52 double sum, norm_tmp = 0.0, norm_b = 0.0;
53 double scaled_res_norm, res_norm, *tmp, max_norm = 0.0;
54
55
56/* Computes the residual
57
58 res = || b - A*x ||
59
60 where x and b are vectors and A is a sparse matrix stored
61 in MSR format. */
62
63/* --------------------------
64 First executable statement
65 -------------------------- */
66
67 /* Create tmp workspace */
68 tmp = (double *) calloc(m,sizeof(double));
69
70/* .....initialize soln */
71
72 for (i = 0; i < m; i++)
73 tmp[i] = b[i] - val[i] * xlocal[i];
74
75/* .....do a series of SPDOTs (sparse dot products) */
76
77 for (i = 0; i <m ; i++)
78 {
79 jbgn = indx[i];
80 jend = indx[i + 1];
81 sum = 0.0;
82
83 for (j = jbgn; j < jend; j++)
84 sum += val[j] * x[indx[j]];
85
86 tmp[i] -= sum;
87 max_norm = max(fabs(tmp[i]),max_norm);
88 norm_tmp += tmp[i]*tmp[i];
89 norm_b += b[i]*b[i];
90 }
91
92 res_norm = sqrt(norm_tmp);
93 printf("\n\nMax norm of residual = %12.4g\n",max_norm);
94 printf( "Two norm of residual = %12.4g\n",res_norm);
95 if (norm_b > 1.0E-7)
96 {
97 scaled_res_norm = res_norm/sqrt(norm_b);
98 printf( "Scaled two norm of residual = %12.4g\n",scaled_res_norm);
99 }
100 free((void *) tmp);
101
102 return(scaled_res_norm);
103
104} /* smsrres */
105
double smsrres(int m, int n, double *val, int *indx, double *xlocal, double *x, double *b)
Definition smsrres.c:47
#define max(x, y)
Definition smsrres.c:46