Ifpack Package Browser (Single Doxygen Collection) Development
svbrres.c
Go to the documentation of this file.
00001 /*@HEADER
00002 // ***********************************************************************
00003 // 
00004 //       Ifpack: Object-Oriented Algebraic Preconditioner Package
00005 //                 Copyright (2002) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 //@HEADER
00028 */
00029 
00030 #include <stdlib.h>
00031 #include <stdio.h>
00032 #include <math.h>
00033 #include "spblas.h"
00034 #define max(x,y) (( x > y ) ? x : y)     /* max function  */
00035 
00036 double svbrres (int m, int n, int m_blk,
00037     double *val, int *indx, int *bindx, int *rpntr,
00038     int *cpntr, int *bpntrb, int *bpntre,
00039     double *x, double *b)
00040 {
00041     int i, j, jbgn, jend, ione = 1;
00042     double sum, norm_tmp = 0.0, norm_b = 0.0;
00043     double scaled_res_norm, res_norm, *tmp, max_norm = 0.0;
00044     SPBLASMAT  A;
00045 
00046 
00047 
00048 /*     Computes the residual
00049 
00050                       res = || b - A*x ||
00051 
00052        where x and b are vectors and A is a sparse matrix stored
00053        in MSR format. */
00054 
00055 /*     -------------------------- 
00056        First executable statement 
00057        -------------------------- */
00058     /* Create sparse matrix handle */
00059     cblas_duscr_vbr(m_blk, val, indx, bindx, rpntr, cpntr, bpntrb, bpntre, &A);
00060     
00061 
00062     /* Create tmp workspace, set to b */
00063 
00064     tmp = (double *) calloc(m,sizeof(double));
00065 
00066     for (i = 0; i < m; i++) tmp[i] = b[i];
00067 
00068     /* Call DUSMM to compute residual (in tmp) */
00069 
00070     cblas_dusmm(m_blk, 1, n, -1.0, &A, x, m, 1.0, tmp, m);
00071 
00072     for (i = 0; i <m ; i++) 
00073       {
00074   max_norm = max(fabs(tmp[i]),max_norm);
00075   norm_tmp += tmp[i]*tmp[i];
00076   norm_b += b[i]*b[i];
00077       }
00078    
00079     res_norm = sqrt(norm_tmp);
00080     scaled_res_norm = res_norm/sqrt(norm_b);
00081     printf("\n\nMax norm of residual        = %12.4g\n",max_norm);
00082     printf(    "Two norm of residual        = %12.4g\n",res_norm);
00083     if (norm_b > 1.0E-7) 
00084       {
00085   scaled_res_norm = res_norm/sqrt(norm_b);
00086   printf(    "Scaled two norm of residual = %12.4g\n",scaled_res_norm);
00087       }
00088     /* Compute residual statistics */
00089     /*      if (res_norm > 0.2 )
00090     cblas_dusmm_dump("/u1/mheroux/dump_file",
00091          m_blk, 1, n, -1.0, &A, x, n, 1.0, b, m);
00092    for (i=0; i<m_blk; i++)
00093       {
00094   printf("***** Row %d *******\n",i);
00095   printf("bpntrb[%d] = %d\n",i,bpntrb[i]);
00096   printf("bpntre[%d] = %d\n",i,bpntre[i]);
00097   printf("rpntr[%d] = %d\n",i,rpntr[i]);
00098   for (j=bpntrb[i]; j<bpntre[i]; j++)
00099     {
00100       printf("bindx[%d] = %d\n",j,bindx[j]);
00101       printf("indx[%d] = %d\n",j,indx[j]);
00102     }
00103   
00104     
00105       }
00106   printf("rpntr[%d] = %d\n",m_blk,rpntr[m_blk]);
00107   j = bpntre[m_blk-1];
00108   printf("bindx[%d] = %d\n",j,bindx[j]);
00109   printf("indx[%d] = %d\n",j,indx[j]);
00110     printf("val[indx[bpntrb[m_blk-1]]  ] = %12.4g\n",val[indx[bpntrb[m_blk-1]]  ]);
00111     printf("val[indx[bpntrb[m_blk-1]]+1] = %12.4g\n",val[indx[bpntrb[m_blk-1]]+1]);
00112     printf("val[indx[bpntrb[m_blk-1]]+2] = %12.4g\n",val[indx[bpntrb[m_blk-1]]+2]);
00113 
00114     for (i = 0; i <m ; i++) 
00115       {
00116   printf("tmp[%d] = %12.4g\n",i,tmp[i]);
00117   printf("  x[%d] = %12.4g\n",i,  x[i]);
00118   printf("  b[%d] = %12.4g\n",i,  b[i]);
00119       }
00120     */
00121     free((void *) tmp);
00122 
00123     return(res_norm);
00124 
00125 } /* svbrres */
00126 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines