Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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
49double 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