read_hb.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 "iohb.h"
00033 #include "prototypes.h"
00034 
00035 void read_hb(char *data_file,
00036         int *N_global, int *n_nonzeros, 
00037         double **val, int **bindx,
00038         double **x, double **b, double **bt, double **xexact)
00039 #undef DEBUG 
00040 
00041 {
00042   FILE *in_file ;
00043   char Title[73], Key[9], Rhstype[4];
00044   char Type[4] = "XXX\0";
00045   char Ptrfmt[17], Indfmt[17], Valfmt[21], Rhsfmt[21];
00046   int Ptrcrd, Indcrd, Valcrd, Rhscrd;
00047 
00048   int i, n_entries, N_columns, Nrhs;
00049   int ii, jj ;
00050   int kk = 0;
00051   int isym;
00052   int ione = 1;
00053   double value, res;
00054   double *cnt;
00055   int *pntr, *indx1, *pntr1;
00056   double *val1;
00057 
00058   
00059       in_file = fopen( data_file, "r");
00060       if (in_file == NULL)
00061   {
00062     printf("Error: Cannot open file: %s\n",data_file);
00063     exit(1);
00064   }
00065 
00066       /* Get information about the array stored in the file specified in the  */
00067       /* argument list:                                                       */
00068 
00069       printf("Reading matrix info from %s...\n",data_file);
00070       
00071       in_file = fopen( data_file, "r");
00072       if (in_file == NULL)
00073   {
00074     printf("Error: Cannot open file: %s\n",data_file);
00075     exit(1);
00076   }
00077 
00078       readHB_header(in_file, Title, Key, Type, N_global, &N_columns, 
00079         &n_entries, &Nrhs,
00080         Ptrfmt, Indfmt, Valfmt, Rhsfmt,
00081         &Ptrcrd, &Indcrd, &Valcrd, &Rhscrd, Rhstype);
00082       fclose(in_file);
00083 
00084       if (Nrhs < 0 ) Nrhs = 0;
00085 
00086       printf("***************************************************************\n");
00087       printf("Matrix in file %s is %d x %d, \n",data_file, *N_global, N_columns);
00088       printf("with %d nonzeros with type %3s;\n", n_entries, Type);
00089       printf("***************************************************************\n");
00090       printf("Title: %72s\n",Title);
00091       printf("***************************************************************\n");
00092       /*Nrhs = 0; */
00093       printf("%d right-hand-side(s) available.\n",Nrhs);
00094 
00095       if (Type[0] != 'R') perror("Can only handle real valued matrices");
00096       isym = 0;
00097       if (Type[1] == 'S') 
00098   {
00099     printf("Converting symmetric matrix to nonsymmetric storage\n");
00100     n_entries = 2*n_entries - N_columns;
00101     isym = 1;
00102   }
00103       if (Type[2] != 'A') perror("Can only handle assembled matrices");
00104       if (N_columns != *N_global) perror("Matrix dimensions must be the same");
00105       *n_nonzeros = n_entries;
00106       
00107       /* Read the matrix information, generating the associated storage arrays  */
00108       printf("Reading the matrix from %s...\n",data_file);
00109 
00110       /* Allocate space.  Note that we add extra storage in case of zero
00111    diagonals.  This is necessary for conversion to MSR format. */
00112 
00113       pntr   = (int    *) calloc(N_columns+1,sizeof(int)) ;
00114       *bindx = (int    *) calloc(n_entries+N_columns+1,sizeof(int)) ;
00115       *val   = (double *) calloc(n_entries+N_columns+1,sizeof(double)) ;
00116       
00117       readHB_mat_double(data_file, pntr, *bindx, *val);
00118 
00119       /* If a rhs is specified in the file, read one, 
00120    generating the associate storage */
00121       if (Nrhs > 0 && Rhstype[2] =='X')
00122   {
00123     printf("Reading right-hand-side vector(s) from %s...\n",data_file);
00124     *b = (double *) calloc(N_columns,sizeof(double));
00125     readHB_aux_double(data_file, 'F', (*b));
00126     printf("Reading exact solution  vector(s) from %s...\n",data_file);
00127     *xexact = (double *) calloc(N_columns,sizeof(double));
00128         readHB_aux_double(data_file, 'X', (*xexact));
00129 
00130   }
00131       else
00132   {
00133     
00134     /* Set Xexact to a random vector */
00135 
00136     printf("Setting  random exact solution  vector\n");
00137     *xexact = (double *) calloc(N_columns,sizeof(double));
00138     
00139     for (i=0;i<*N_global;i++)  (*xexact)[i] = drand48();
00140     
00141     /* Compute b to match xexact */
00142     
00143     *b = (double   *) calloc(N_columns,sizeof(double)) ;
00144     if ((*b) == NULL) perror("Error: Not enough space to create rhs");
00145  
00146     scscmv (isym, 1, N_columns, N_columns, (*val), (*bindx), pntr, (*xexact), (*b));
00147 
00148   }
00149 
00150       /* Compute residual using CSC format */
00151 
00152       for (i = 0; i <= *N_global; i++) pntr[i]--;
00153       for (i = 0; i <= n_entries; i++) (*bindx)[i]--;
00154       res = scscres(isym, *N_global, *N_global, (*val), (*bindx), pntr, 
00155         (*xexact), (*b));
00156       printf(
00157         "The residual using CSC format and exact solution is %12.4g\n",
00158         res);
00159       for (i = 0; i <= *N_global; i++) pntr[i]++;
00160       for (i = 0; i <= n_entries; i++) (*bindx)[i]++;
00161 
00162       
00163       /* Set initial guess to zero */
00164       
00165       *x = (double   *) calloc((*N_global),sizeof(double)) ;
00166       
00167       if ((*x) == NULL) 
00168   perror("Error: Not enough space to create guess");
00169       
00170       
00171       /* Set RHS to a random vector, initial guess to zero */
00172       for (i=0;i<*N_global;i++) (*x)[i] = 0.0;
00173       
00174       
00175       /* Allocate temporary space */
00176       
00177       pntr1 = (int   *) calloc(N_columns+1,sizeof(int)) ;
00178       indx1 = (int   *) calloc(n_entries+N_columns+1,sizeof(int)) ;
00179       val1 = (double *) calloc(n_entries+N_columns+1,sizeof(double)) ;
00180       
00181       
00182       /* Convert in the following way:
00183    - CSC to CSR 
00184    - CSR to MSR
00185       */
00186       csrcsc_(N_global,&ione,&ione,
00187         *val,*bindx,pntr,
00188         val1,indx1,pntr1);
00189 
00190       /* Compute bt = A(trans)*xexact */
00191 
00192       *bt = (double   *) calloc(N_columns,sizeof(double)) ;
00193       
00194       scscmv (isym, 1, N_columns, N_columns, val1, indx1, pntr1, (*xexact), (*bt));
00195 
00196       if (Type[1] == 'S') 
00197   {
00198     int *indu;
00199     int ierr;
00200     indu = indx1+n_entries; /* Use end of bindx for workspace */
00201     ssrcsr_(&N_columns, val1, indx1, pntr1, &n_entries,
00202           val1, indx1, pntr1, indu, &ierr);
00203     if (ierr !=0 ) 
00204       {
00205       printf(" Error in converting from symmetric form\n  IERR = %d\n",ierr);
00206       abort();
00207       }
00208   }
00209 
00210       csrmsr_(N_global,val1,indx1,pntr1,
00211         *val,*bindx,
00212         *val,*bindx);
00213       
00214       /* Recompute number of nonzeros in case there were zero diagonals */
00215       
00216       *n_nonzeros = (*bindx)[*N_global] - 2; /* Still in Fortran mode so -2 */
00217       
00218       /* Finally, convert bindx vectors to zero base */
00219       
00220       for (i=0;i<*n_nonzeros+1;i++) (*bindx)[i] -= 1;
00221       
00222 
00223       printf("The residual using MSR format and exact solution is %12.4g\n",
00224         smsrres (*N_global, *N_global, (*val), (*bindx), 
00225            (*xexact), (*xexact), (*b)));
00226 
00227       /* Release unneeded space */
00228 
00229       free((void *) val1);
00230       free((void *) indx1);
00231       free((void *) pntr1);
00232       free((void *) pntr);
00233   
00234   /* end read_hb */
00235 }

Generated on Thu Sep 18 12:37:22 2008 for Ifpack Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1