Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
scscres.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 scscres (int isym, int m, int n,
48 double *val, int *indx, int *pntr,
49 double *x, double *b)
50{
51 int i, j, ibgn, iend, ione = 1;
52 double 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];
74
75/* .....do a series of SPAXPYs (sparse saxpys) */
76
77 for (j = 0; j < n ; j++)
78 {
79 ibgn = pntr[j];
80 iend = pntr[j + 1];
81
82 for (i = ibgn; i < iend; i++)
83 {
84 tmp[indx[i]] -= val[i] * x[j];
85 if (indx[i] != j && isym) tmp[j] -= val[i]*x[indx[i]];
86 }
87 }
88 for (i = 0; i < m; i++)
89 {
90 max_norm = max(fabs(tmp[i]),max_norm);
91 norm_tmp += tmp[i]*tmp[i];
92 norm_b += b[i]*b[i];
93 }
94
95 res_norm = sqrt(norm_tmp);
96 printf("\n\nMax norm of residual = %12.4g\n",max_norm);
97 printf( "Two norm of residual = %12.4g\n",res_norm);
98 if (norm_b > 1.0E-7)
99 {
100 scaled_res_norm = res_norm/sqrt(norm_b);
101 printf( "Scaled two norm of residual = %12.4g\n",scaled_res_norm);
102 }
103
104 free((void *) tmp);
105
106 return(scaled_res_norm);
107
108} /* scscres */
109
double scscres(int isym, int m, int n, double *val, int *indx, int *pntr, double *x, double *b)
Definition: scscres.c:47
#define max(x, y)
Definition: scscres.c:46