49 int *N_global,
int *n_nonzeros,
50 double **val,
int **bindx,
51 double **x,
double **b,
double **bt,
double **xexact)
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;
61 int i, n_entries, N_columns, Nrhs;
68 int *pntr, *indx1, *pntr1;
72 in_file = fopen( data_file,
"r");
75 printf(
"Error: Cannot open file: %s\n",data_file);
82 printf(
"Reading matrix info from %s...\n",data_file);
84 in_file = fopen( data_file,
"r");
87 printf(
"Error: Cannot open file: %s\n",data_file);
91 readHB_header(in_file, Title, Key, Type, N_global, &N_columns,
93 Ptrfmt, Indfmt, Valfmt, Rhsfmt,
94 &Ptrcrd, &Indcrd, &Valcrd, &Rhscrd, Rhstype);
97 if (Nrhs < 0 ) Nrhs = 0;
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");
106 printf(
"%d right-hand-side(s) available.\n",Nrhs);
108 if (Type[0] !=
'R')
perror(
"Can only handle real valued matrices");
112 printf(
"Converting symmetric matrix to nonsymmetric storage\n");
113 n_entries = 2*n_entries - N_columns;
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;
121 printf(
"Reading the matrix from %s...\n",data_file);
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)) ;
134 if (Nrhs > 0 && Rhstype[2] ==
'X')
136 printf(
"Reading right-hand-side vector(s) from %s...\n",data_file);
137 *b = (
double *) calloc(N_columns,
sizeof(
double));
139 printf(
"Reading exact solution vector(s) from %s...\n",data_file);
140 *xexact = (
double *) calloc(N_columns,
sizeof(
double));
149 printf(
"Setting random exact solution vector\n");
150 *xexact = (
double *) calloc(N_columns,
sizeof(
double));
152 for (i=0;i<*N_global;i++) (*xexact)[i] = drand48();
156 *b = (
double *) calloc(N_columns,
sizeof(
double)) ;
157 if ((*b) == NULL)
perror(
"Error: Not enough space to create rhs");
159 scscmv (isym, 1, N_columns, N_columns, (*val), (*bindx), pntr, (*xexact), (*b));
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,
170 "The residual using CSC format and exact solution is %12.4g\n",
172 for (i = 0; i <= *N_global; i++) pntr[i]++;
173 for (i = 0; i <= n_entries; i++) (*bindx)[i]++;
178 *x = (
double *) calloc((*N_global),
sizeof(double)) ;
181 perror(
"Error: Not enough space to create guess");
185 for (i=0;i<*N_global;i++) (*x)[i] = 0.0;
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)) ;
199 csrcsc_(N_global,&ione,&ione,
205 *bt = (
double *) calloc(N_columns,
sizeof(
double)) ;
207 scscmv (isym, 1, N_columns, N_columns, val1, indx1, pntr1, (*xexact), (*bt));
213 indu = indx1+n_entries;
214 ssrcsr_(&N_columns, val1, indx1, pntr1, &n_entries,
215 val1, indx1, pntr1, indu, &ierr);
218 printf(
" Error in converting from symmetric form\n IERR = %d\n",ierr);
223 csrmsr_(N_global,val1,indx1,pntr1,
229 *n_nonzeros = (*bindx)[*N_global] - 2;
233 for (i=0;i<*n_nonzeros+1;i++) (*bindx)[i] -= 1;
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)));
243 free((
void *) indx1);
244 free((
void *) pntr1);
double scscres(int isym, int m, int n, double *val, int *indx, int *pntr, double *x, double *b)
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)
double smsrres(int m, int n, double *val, int *indx, double *xlocal, double *x, double *b)
int readHB_mat_double(const char *filename, int colptr[], int rowind[], double val[])
void read_hb(char *data_file, int *N_global, int *n_nonzeros, double **val, int **bindx, double **x, double **b, double **bt, double **xexact)
void scscmv(int isym, int indexbase, int m, int n, double *val, int *indx, int *pntr, double *x, double *b)
int readHB_aux_double(const char *filename, const char AuxType, double b[])