Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
read_hb.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 "iohb.h"
46#include "prototypes.h"
47
48void read_hb(char *data_file,
49 int *N_global, int *n_nonzeros,
50 double **val, int **bindx,
51 double **x, double **b, double **bt, double **xexact)
52#undef DEBUG
53
54{
55 FILE *in_file ;
56 char Title[73], Key[9], Rhstype[4];
57 char Type[4] = "XXX\0";
58 char Ptrfmt[17], Indfmt[17], Valfmt[21], Rhsfmt[21];
59 int Ptrcrd, Indcrd, Valcrd, Rhscrd;
60
61 int i, n_entries, N_columns, Nrhs;
62 int ii, jj ;
63 int kk = 0;
64 int isym;
65 int ione = 1;
66 double value, res;
67 double *cnt;
68 int *pntr, *indx1, *pntr1;
69 double *val1;
70
71
72 in_file = fopen( data_file, "r");
73 if (in_file == NULL)
74 {
75 printf("Error: Cannot open file: %s\n",data_file);
76 exit(1);
77 }
78
79 /* Get information about the array stored in the file specified in the */
80 /* argument list: */
81
82 printf("Reading matrix info from %s...\n",data_file);
83
84 in_file = fopen( data_file, "r");
85 if (in_file == NULL)
86 {
87 printf("Error: Cannot open file: %s\n",data_file);
88 exit(1);
89 }
90
91 readHB_header(in_file, Title, Key, Type, N_global, &N_columns,
92 &n_entries, &Nrhs,
93 Ptrfmt, Indfmt, Valfmt, Rhsfmt,
94 &Ptrcrd, &Indcrd, &Valcrd, &Rhscrd, Rhstype);
95 fclose(in_file);
96
97 if (Nrhs < 0 ) Nrhs = 0;
98
99 printf("***************************************************************\n");
100 printf("Matrix in file %s is %d x %d, \n",data_file, *N_global, N_columns);
101 printf("with %d nonzeros with type %3s;\n", n_entries, Type);
102 printf("***************************************************************\n");
103 printf("Title: %72s\n",Title);
104 printf("***************************************************************\n");
105 /*Nrhs = 0; */
106 printf("%d right-hand-side(s) available.\n",Nrhs);
107
108 if (Type[0] != 'R') perror("Can only handle real valued matrices");
109 isym = 0;
110 if (Type[1] == 'S')
111 {
112 printf("Converting symmetric matrix to nonsymmetric storage\n");
113 n_entries = 2*n_entries - N_columns;
114 isym = 1;
115 }
116 if (Type[2] != 'A') perror("Can only handle assembled matrices");
117 if (N_columns != *N_global) perror("Matrix dimensions must be the same");
118 *n_nonzeros = n_entries;
119
120 /* Read the matrix information, generating the associated storage arrays */
121 printf("Reading the matrix from %s...\n",data_file);
122
123 /* Allocate space. Note that we add extra storage in case of zero
124 diagonals. This is necessary for conversion to MSR format. */
125
126 pntr = (int *) calloc(N_columns+1,sizeof(int)) ;
127 *bindx = (int *) calloc(n_entries+N_columns+1,sizeof(int)) ;
128 *val = (double *) calloc(n_entries+N_columns+1,sizeof(double)) ;
129
130 readHB_mat_double(data_file, pntr, *bindx, *val);
131
132 /* If a rhs is specified in the file, read one,
133 generating the associate storage */
134 if (Nrhs > 0 && Rhstype[2] =='X')
135 {
136 printf("Reading right-hand-side vector(s) from %s...\n",data_file);
137 *b = (double *) calloc(N_columns,sizeof(double));
138 readHB_aux_double(data_file, 'F', (*b));
139 printf("Reading exact solution vector(s) from %s...\n",data_file);
140 *xexact = (double *) calloc(N_columns,sizeof(double));
141 readHB_aux_double(data_file, 'X', (*xexact));
142
143 }
144 else
145 {
146
147 /* Set Xexact to a random vector */
148
149 printf("Setting random exact solution vector\n");
150 *xexact = (double *) calloc(N_columns,sizeof(double));
151
152 for (i=0;i<*N_global;i++) (*xexact)[i] = drand48();
153
154 /* Compute b to match xexact */
155
156 *b = (double *) calloc(N_columns,sizeof(double)) ;
157 if ((*b) == NULL) perror("Error: Not enough space to create rhs");
158
159 scscmv (isym, 1, N_columns, N_columns, (*val), (*bindx), pntr, (*xexact), (*b));
160
161 }
162
163 /* Compute residual using CSC format */
164
165 for (i = 0; i <= *N_global; i++) pntr[i]--;
166 for (i = 0; i <= n_entries; i++) (*bindx)[i]--;
167 res = scscres(isym, *N_global, *N_global, (*val), (*bindx), pntr,
168 (*xexact), (*b));
169 printf(
170 "The residual using CSC format and exact solution is %12.4g\n",
171 res);
172 for (i = 0; i <= *N_global; i++) pntr[i]++;
173 for (i = 0; i <= n_entries; i++) (*bindx)[i]++;
174
175
176 /* Set initial guess to zero */
177
178 *x = (double *) calloc((*N_global),sizeof(double)) ;
179
180 if ((*x) == NULL)
181 perror("Error: Not enough space to create guess");
182
183
184 /* Set RHS to a random vector, initial guess to zero */
185 for (i=0;i<*N_global;i++) (*x)[i] = 0.0;
186
187
188 /* Allocate temporary space */
189
190 pntr1 = (int *) calloc(N_columns+1,sizeof(int)) ;
191 indx1 = (int *) calloc(n_entries+N_columns+1,sizeof(int)) ;
192 val1 = (double *) calloc(n_entries+N_columns+1,sizeof(double)) ;
193
194
195 /* Convert in the following way:
196 - CSC to CSR
197 - CSR to MSR
198 */
199 csrcsc_(N_global,&ione,&ione,
200 *val,*bindx,pntr,
201 val1,indx1,pntr1);
202
203 /* Compute bt = A(trans)*xexact */
204
205 *bt = (double *) calloc(N_columns,sizeof(double)) ;
206
207 scscmv (isym, 1, N_columns, N_columns, val1, indx1, pntr1, (*xexact), (*bt));
208
209 if (Type[1] == 'S')
210 {
211 int *indu;
212 int ierr;
213 indu = indx1+n_entries; /* Use end of bindx for workspace */
214 ssrcsr_(&N_columns, val1, indx1, pntr1, &n_entries,
215 val1, indx1, pntr1, indu, &ierr);
216 if (ierr !=0 )
217 {
218 printf(" Error in converting from symmetric form\n IERR = %d\n",ierr);
219 abort();
220 }
221 }
222
223 csrmsr_(N_global,val1,indx1,pntr1,
224 *val,*bindx,
225 *val,*bindx);
226
227 /* Recompute number of nonzeros in case there were zero diagonals */
228
229 *n_nonzeros = (*bindx)[*N_global] - 2; /* Still in Fortran mode so -2 */
230
231 /* Finally, convert bindx vectors to zero base */
232
233 for (i=0;i<*n_nonzeros+1;i++) (*bindx)[i] -= 1;
234
235
236 printf("The residual using MSR format and exact solution is %12.4g\n",
237 smsrres (*N_global, *N_global, (*val), (*bindx),
238 (*xexact), (*xexact), (*b)));
239
240 /* Release unneeded space */
241
242 free((void *) val1);
243 free((void *) indx1);
244 free((void *) pntr1);
245 free((void *) pntr);
246
247 /* end read_hb */
248}
#define perror(str)
Definition cc_main.cc:55
int readHB_header(FILE *in_file, char *Title, char *Key, char *Type, int *Nrow, int *Ncol, int *Nnzero, int *Nrhs, char *Ptrfmt, char *Indfmt, char *Valfmt, char *Rhsfmt, int *Ptrcrd, int *Indcrd, int *Valcrd, int *Rhscrd, char *Rhstype)
Definition iohb.c:293
int readHB_aux_double(const char *filename, const char AuxType, double b[])
Definition iohb.c:543
int readHB_mat_double(const char *filename, int colptr[], int rowind[], double val[])
Definition iohb.c:366
void scscmv(int isym, int indexbase, int m, int n, double *val, int *indx, int *pntr, double *x, double *b)
Definition scscmv.c:43
double scscres(int isym, int m, int n, double *val, int *indx, int *pntr, double *x, double *b)
Definition scscres.c:47
double smsrres(int m, int n, double *val, int *indx, double *xlocal, double *x, double *b)
Definition smsrres.c:47
void read_hb(char *data_file, int *N_global, int *n_nonzeros, double **val, int **bindx, double **x, double **b, double **bt, double **xexact)
Definition read_hb.c:48