Ifpack Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
svbrres.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 #include "spblas.h"
47 #define max(x,y) (( x > y ) ? x : y) /* max function */
48 
49 double svbrres (int m, int n, int m_blk,
50  double *val, int *indx, int *bindx, int *rpntr,
51  int *cpntr, int *bpntrb, int *bpntre,
52  double *x, double *b)
53 {
54  int i, j, jbgn, jend, ione = 1;
55  double sum, norm_tmp = 0.0, norm_b = 0.0;
56  double scaled_res_norm, res_norm, *tmp, max_norm = 0.0;
57  SPBLASMAT A;
58 
59 
60 
61 /* Computes the residual
62 
63  res = || b - A*x ||
64 
65  where x and b are vectors and A is a sparse matrix stored
66  in MSR format. */
67 
68 /* --------------------------
69  First executable statement
70  -------------------------- */
71  /* Create sparse matrix handle */
72  cblas_duscr_vbr(m_blk, val, indx, bindx, rpntr, cpntr, bpntrb, bpntre, &A);
73 
74 
75  /* Create tmp workspace, set to b */
76 
77  tmp = (double *) calloc(m,sizeof(double));
78 
79  for (i = 0; i < m; i++) tmp[i] = b[i];
80 
81  /* Call DUSMM to compute residual (in tmp) */
82 
83  cblas_dusmm(m_blk, 1, n, -1.0, &A, x, m, 1.0, tmp, m);
84 
85  for (i = 0; i <m ; i++)
86  {
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  scaled_res_norm = res_norm/sqrt(norm_b);
94  printf("\n\nMax norm of residual = %12.4g\n",max_norm);
95  printf( "Two norm of residual = %12.4g\n",res_norm);
96  if (norm_b > 1.0E-7)
97  {
98  scaled_res_norm = res_norm/sqrt(norm_b);
99  printf( "Scaled two norm of residual = %12.4g\n",scaled_res_norm);
100  }
101  /* Compute residual statistics */
102  /* if (res_norm > 0.2 )
103  cblas_dusmm_dump("/u1/mheroux/dump_file",
104  m_blk, 1, n, -1.0, &A, x, n, 1.0, b, m);
105  for (i=0; i<m_blk; i++)
106  {
107  printf("***** Row %d *******\n",i);
108  printf("bpntrb[%d] = %d\n",i,bpntrb[i]);
109  printf("bpntre[%d] = %d\n",i,bpntre[i]);
110  printf("rpntr[%d] = %d\n",i,rpntr[i]);
111  for (j=bpntrb[i]; j<bpntre[i]; j++)
112  {
113  printf("bindx[%d] = %d\n",j,bindx[j]);
114  printf("indx[%d] = %d\n",j,indx[j]);
115  }
116 
117 
118  }
119  printf("rpntr[%d] = %d\n",m_blk,rpntr[m_blk]);
120  j = bpntre[m_blk-1];
121  printf("bindx[%d] = %d\n",j,bindx[j]);
122  printf("indx[%d] = %d\n",j,indx[j]);
123  printf("val[indx[bpntrb[m_blk-1]] ] = %12.4g\n",val[indx[bpntrb[m_blk-1]] ]);
124  printf("val[indx[bpntrb[m_blk-1]]+1] = %12.4g\n",val[indx[bpntrb[m_blk-1]]+1]);
125  printf("val[indx[bpntrb[m_blk-1]]+2] = %12.4g\n",val[indx[bpntrb[m_blk-1]]+2]);
126 
127  for (i = 0; i <m ; i++)
128  {
129  printf("tmp[%d] = %12.4g\n",i,tmp[i]);
130  printf(" x[%d] = %12.4g\n",i, x[i]);
131  printf(" b[%d] = %12.4g\n",i, b[i]);
132  }
133  */
134  free((void *) tmp);
135 
136  return(res_norm);
137 
138 } /* svbrres */
139 
double svbrres(int m, int n, int m_blk, double *val, int *indx, int *bindx, int *rpntr, int *cpntr, int *bpntrb, int *bpntre, double *x, double *b)
Definition: svbrres.c:49
#define max(x, y)
Definition: svbrres.c:47