Dear PETSc Team,
Sorry to bother! My name is Hemeng Wang, and I am currently learning the use of PETSc software package. I am confused while calculating the norm of residue. I calculated residue norm by myself with: ``` PetscCall(VecNorm(b, NORM_2, &norm_b)); // (main.c, line 74) PetscCall(VecDuplicate(b, &u)); // (main.c, line 105) PetscCall(MatMult(A, x, u)); PetscCall(VecAXPY(b, -1.0, u)); PetscCall(VecNorm(b, NORM_2, &norm_delta)); ``` and check the (norm_delta) / (norm_b). It seems not the `atol` which set by `KSPSetTolerances()`. (I set atol as 1e-12, but got 5e-7 finally) (options: -ksp_type cg -pc_type gamg -ksp_converged_reason -ksp_norm_type unpreconditioned -ksp_monitor_true_residual) I also check the soure code of `KSPSolve_CG` in `petsc/src/ksp/ksp/impls/cg/cg.c`. And could not figure out where is the difference. I will really really appreciated if someone can explain the calculation of `resid norm` in petsc. And where is my mistake. To provide you with more context, here are the source code about my implementation. And the output of my test. main.c Main code of my program mmio.h mmloader.h Headers for matrix read Makefile For compiling, same as sample provided task.sh A script for running program in `slurm` slurm-4803840.out Output of my test Thank you very much for your time and attention. I greatly appreciate your support and look forward to hearing from you soon. Best regards, Hemeng Wang
static char help[] = "Solves a linear system in parallel with KSP.\n\
Input parameters include:\n\
-view_exact_sol : write exact solution vector to stdout\n\
-m <mesh_x> : number of mesh points in x-direction\n\
-n <mesh_y> : number of mesh points in y-direction\n\n";
/*
Include "petscksp.h" so that we can use KSP solvers.
*/
#include <petscksp.h>
#include <sys/time.h>
#include <float.h>
#include "mmloader.h"
int main(int argc, char **args)
{
Vec x, b; /* approx solution, RHS */
Vec u;
Mat A; /* linear system matrix */
KSP ksp; /* linear solver context */
PetscReal residue; /* residue of solution error */
PetscInt m, n, its;
PetscReal norm_delta, norm_b; /* norm of solution error */
char filename[PETSC_MAX_PATH_LEN], bfilename[PETSC_MAX_PATH_LEN];
PetscFunctionBeginUser;
PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
PetscCall(PetscOptionsGetString(NULL, NULL, "-MTX", filename,
PETSC_MAX_PATH_LEN, NULL));
PetscCall(PetscOptionsGetString(NULL, NULL, "-RHS", bfilename,
PETSC_MAX_PATH_LEN, NULL));
PetscCall(PetscPrintf(PETSC_COMM_WORLD, "mat name: %s\n", filename));
PetscCall(PetscPrintf(PETSC_COMM_WORLD, "rhs name: %s\n", bfilename));
// PetscCall(MatCreateFromMTX(&A, filename, PETSC_TRUE));
PetscCall(load_binary_mat(&A, filename));
PetscCall(MatGetSize(A, &m, &n));
// PetscCall(VecCreateFromRHS(&b, bfilename));
PetscCall(load_binary_rhs(&b, bfilename));
PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
PetscCall(VecSetSizes(x, PETSC_DECIDE, n));
PetscCall(VecSetFromOptions(x));
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Create the linear solver and set various options
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
/*
Set operators. Here the matrix that defines the linear system
also serves as the preconditioning matrix.
*/
PetscCall(KSPSetOperators(ksp, A, A));
/*
Set linear solver defaults for this problem (optional).
- By extracting the KSP and PC contexts from the KSP context,
we can then directly call any KSP and PC routines to set
various options.
- The following two statements are optional; all of these
parameters could alternatively be specified at runtime via
KSPSetFromOptions(). All of these defaults can be
overridden at runtime, as indicated below.
*/
PetscCall(VecNorm(b, NORM_2, &norm_b));
// PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Atol: %lf\n", 1.e-10 * norm_b));
PetscCall(KSPSetTolerances(ksp, 1e-12, DBL_MIN, PETSC_DEFAULT,
PETSC_DEFAULT));
/*
Set runtime options, e.g.,
-ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
These options will override those specified above as long as
KSPSetFromOptions() is called _after_ any other customization
routines.
*/
PetscCall(KSPSetFromOptions(ksp));
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Solve the linear system
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
struct timeval t_start, t_stop;
gettimeofday(&t_start, NULL);
PetscCall(KSPSolve(ksp, b, x));
gettimeofday(&t_stop, NULL);
double total_time = (t_stop.tv_sec - t_start.tv_sec) * 1000.0 +
(t_stop.tv_usec - t_start.tv_usec) / 1000.0;
PetscCall(PetscPrintf(PETSC_COMM_WORLD, "KSPSolve Time: %lf ms\n",
total_time));
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Check the solution and clean up
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
PetscCall(VecDuplicate(b, &u));
PetscCall(MatMult(A, x, u));
PetscCall(VecAXPY(b, -1.0, u));
PetscCall(VecNorm(b, NORM_2, &norm_delta));
PetscCall(KSPGetResidualNorm(ksp, &residue));
PetscCall(KSPGetIterationNumber(ksp, &its));
PetscCall(PetscPrintf(PETSC_COMM_WORLD, "norm_delta %lf, norm_b %lf \n",
norm_delta, norm_b));
PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g iterations %"
PetscInt_FMT "\n", (norm_delta) / (norm_b), its));
PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Residue of error %g iterations %"
PetscInt_FMT "\n", (double)residue, its));
/*
Free work space. All PETSc objects should be destroyed when they
are no longer needed.
*/
// PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Free work space\n"));
PetscCall(KSPDestroy(&ksp));
PetscCall(VecDestroy(&x));
PetscCall(VecDestroy(&b));
PetscCall(MatDestroy(&A));
/*
Always call PetscFinalize() before exiting a program. This routine
- finalizes the PETSc libraries as well as MPI
- provides summary and diagnostic information if certain runtime
options are chosen (e.g., -log_view).
*/
PetscCall(PetscFinalize());
return 0;
}
/*TEST
build:
requires: !single
test:
suffix: chebyest_1
args: -m 80 -n 80 -ksp_pc_side right -pc_type ksp -ksp_ksp_type chebyshev
-ksp_ksp_max_it 5 -ksp_ksp_chebyshev_esteig 0.9,0,0,1.1 -ksp_monitor_short
test:
suffix: chebyest_2
args: -m 80 -n 80 -ksp_pc_side right -pc_type ksp -ksp_ksp_type chebyshev
-ksp_ksp_max_it 5 -ksp_ksp_chebyshev_esteig 0.9,0,0,1.1 -ksp_esteig_ksp_type cg
-ksp_monitor_short
test:
args: -ksp_monitor_short -m 5 -n 5 -ksp_gmres_cgs_refinement_type
refine_always
test:
suffix: 2
nsize: 2
args: -ksp_monitor_short -m 5 -n 5 -ksp_gmres_cgs_refinement_type
refine_always
test:
suffix: 3
args: -pc_type sor -pc_sor_symmetric -ksp_monitor_short
-ksp_gmres_cgs_refinement_type refine_always
test:
suffix: 4
args: -pc_type eisenstat -ksp_monitor_short
-ksp_gmres_cgs_refinement_type refine_always
test:
suffix: 5
nsize: 2
args: -ksp_monitor_short -m 5 -n 5 -mat_view draw
-ksp_gmres_cgs_refinement_type refine_always -nox
output_file: output/ex2_2.out
test:
suffix: bjacobi
nsize: 4
args: -pc_type bjacobi -pc_bjacobi_blocks 1 -ksp_monitor_short
-sub_pc_type jacobi -sub_ksp_type gmres
test:
suffix: bjacobi_2
nsize: 4
args: -pc_type bjacobi -pc_bjacobi_blocks 2 -ksp_monitor_short
-sub_pc_type jacobi -sub_ksp_type gmres -ksp_view
test:
suffix: bjacobi_3
nsize: 4
args: -pc_type bjacobi -pc_bjacobi_blocks 4 -ksp_monitor_short
-sub_pc_type jacobi -sub_ksp_type gmres
test:
suffix: qmrcgs
args: -ksp_type qmrcgs -pc_type ilu
output_file: output/ex2_fbcgs.out
test:
suffix: qmrcgs_2
nsize: 3
args: -ksp_type qmrcgs -pc_type bjacobi
output_file: output/ex2_fbcgs_2.out
test:
suffix: fbcgs
args: -ksp_type fbcgs -pc_type ilu
test:
suffix: fbcgs_2
nsize: 3
args: -ksp_type fbcgsr -pc_type bjacobi
test:
suffix: groppcg
args: -ksp_monitor_short -ksp_type groppcg -m 9 -n 9
test:
suffix: mkl_pardiso_cholesky
requires: mkl_pardiso
args: -ksp_type preonly -pc_type cholesky -mat_type sbaij
-pc_factor_mat_solver_type mkl_pardiso
test:
suffix: mkl_pardiso_lu
requires: mkl_pardiso
args: -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mkl_pardiso
test:
suffix: pipebcgs
args: -ksp_monitor_short -ksp_type pipebcgs -m 9 -n 9
test:
suffix: pipecg
args: -ksp_monitor_short -ksp_type pipecg -m 9 -n 9
test:
suffix: pipecgrr
args: -ksp_monitor_short -ksp_type pipecgrr -m 9 -n 9
test:
suffix: pipecr
args: -ksp_monitor_short -ksp_type pipecr -m 9 -n 9
test:
suffix: pipelcg
args: -ksp_monitor_short -ksp_type pipelcg -m 9 -n 9 -pc_type none
-ksp_pipelcg_pipel 2 -ksp_pipelcg_lmax 2
filter: grep -v "sqrt breakdown in iteration"
test:
suffix: sell
args: -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always -m
9 -n 9 -mat_type sell
test:
requires: mumps
suffix: sell_mumps
args: -ksp_type preonly -m 9 -n 12 -mat_type sell -pc_type lu
-pc_factor_mat_solver_type mumps -pc_factor_mat_ordering_type natural
test:
suffix: telescope
nsize: 4
args: -m 100 -n 100 -ksp_converged_reason -pc_type telescope
-pc_telescope_reduction_factor 4 -telescope_pc_type bjacobi
test:
suffix: umfpack
requires: suitesparse
args: -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type umfpack
test:
suffix: spqr
requires: suitesparse
args: -ksp_type preonly -pc_type qr -pc_factor_mat_solver_type spqr
test:
suffix: pc_symmetric
args: -m 10 -n 9 -ksp_converged_reason -ksp_type gmres -ksp_pc_side
symmetric -pc_type cholesky
test:
suffix: pipeprcg
args: -ksp_monitor_short -ksp_type pipeprcg -m 9 -n 9
test:
suffix: pipeprcg_rcw
args: -ksp_monitor_short -ksp_type pipeprcg -recompute_w false -m 9 -n 9
test:
suffix: pipecg2
requires: !defined(PETSC_HAVE_THREADSAFETY)
args: -ksp_monitor_short -ksp_type pipecg2 -m 9 -n 9 -ksp_norm_type
{{preconditioned unpreconditioned natural}}
test:
suffix: pipecg2_2
requires: !defined(PETSC_HAVE_THREADSAFETY)
nsize: 4
args: -ksp_monitor_short -ksp_type pipecg2 -m 15 -n 9 -ksp_norm_type
{{preconditioned unpreconditioned natural}}
test:
suffix: hpddm
nsize: 4
requires: hpddm !__float128 !__fp16
filter: sed -e "s/ iterations 9/ iterations 8/g"
args: -ksp_converged_reason -ksp_type hpddm -ksp_hpddm_precision {{single
double}shared output} -ksp_pc_side {{left right}shared output}
test:
suffix: hpddm___float128
output_file: output/ex2_hpddm.out
nsize: 4
requires: hpddm __float128
filter: sed -e "s/ iterations 9/ iterations 8/g"
args: -ksp_converged_reason -ksp_type hpddm -ksp_hpddm_precision {{double
quadruple}shared output} -ksp_pc_side {{left right}shared output}
test:
suffix: symmetric_pc
nsize: 1
args: -ksp_monitor -ksp_type gmres -pc_type bjacobi -sub_pc_type icc
-ksp_pc_side symmetric
test:
suffix: symmetric_pc2
nsize: 1
args: -ksp_monitor -ksp_type gmres -pc_type bjacobi -sub_pc_type icc
-ksp_pc_side symmetric -pc_bjacobi_blocks 2
TEST*/
Makefile
Description: Binary data
/* Matrix Market I/O library for ANSI C See https://math.nist.gov/MatrixMarket/ for details. */ #ifndef MM_IO_H #define MM_IO_H #include <stdio.h> #define MM_MAX_LINE_LENGTH 1025 #define MatrixMarketBanner "%%MatrixMarket" #define MM_MAX_TOKEN_LENGTH 64 typedef char MM_typecode[4]; char *mm_typecode_to_str(MM_typecode matcode); int mm_read_banner(FILE *f, MM_typecode *matcode); int mm_read_mtx_crd_size(FILE *f, int *M, int *N, int *nz); int mm_read_mtx_array_size(FILE *f, int *M, int *N); int mm_write_banner(FILE *f, MM_typecode matcode); int mm_write_mtx_crd_size(FILE *f, int M, int N, int nz); int mm_write_mtx_array_size(FILE *f, int M, int N); /********************* MM_typecode query functions ***************************/ #define mm_is_matrix(typecode) ((typecode)[0] == 'M') #define mm_is_sparse(typecode) ((typecode)[1] == 'C') #define mm_is_coordinate(typecode) ((typecode)[1] == 'C') #define mm_is_dense(typecode) ((typecode)[1] == 'A') #define mm_is_array(typecode) ((typecode)[1] == 'A') #define mm_is_complex(typecode) ((typecode)[2] == 'C') #define mm_is_real(typecode) ((typecode)[2] == 'R') #define mm_is_pattern(typecode) ((typecode)[2] == 'P') #define mm_is_integer(typecode) ((typecode)[2] == 'I') #define mm_is_symmetric(typecode) ((typecode)[3] == 'S') #define mm_is_general(typecode) ((typecode)[3] == 'G') #define mm_is_skew(typecode) ((typecode)[3] == 'K') #define mm_is_hermitian(typecode) ((typecode)[3] == 'H') int mm_is_valid(MM_typecode matcode); /* too complex for a macro */ /********************* MM_typecode modify functions ***************************/ #define mm_set_matrix(typecode) ((*typecode)[0] = 'M') #define mm_set_coordinate(typecode) ((*typecode)[1] = 'C') #define mm_set_array(typecode) ((*typecode)[1] = 'A') #define mm_set_dense(typecode) mm_set_array(typecode) #define mm_set_sparse(typecode) mm_set_coordinate(typecode) #define mm_set_complex(typecode) ((*typecode)[2] = 'C') #define mm_set_real(typecode) ((*typecode)[2] = 'R') #define mm_set_pattern(typecode) ((*typecode)[2] = 'P') #define mm_set_integer(typecode) ((*typecode)[2] = 'I') #define mm_set_symmetric(typecode) ((*typecode)[3] = 'S') #define mm_set_general(typecode) ((*typecode)[3] = 'G') #define mm_set_skew(typecode) ((*typecode)[3] = 'K') #define mm_set_hermitian(typecode) ((*typecode)[3] = 'H') #define mm_clear_typecode(typecode) ((*typecode)[0] = (*typecode)[1] = (*typecode)[2] = ' ', (*typecode)[3] = 'G') #define mm_initialize_typecode(typecode) mm_clear_typecode(typecode) /********************* Matrix Market error codes ***************************/ #define MM_COULD_NOT_READ_FILE 11 #define MM_PREMATURE_EOF 12 #define MM_NOT_MTX 13 #define MM_NO_HEADER 14 #define MM_UNSUPPORTED_TYPE 15 #define MM_LINE_TOO_LONG 16 #define MM_COULD_NOT_WRITE_FILE 17 /******************** Matrix Market internal definitions ******************** MM_matrix_typecode: 4-character sequence object sparse/ data storage dense type scheme string position: [0] [1] [2] [3] Matrix typecode: M(atrix) C(oord) R(eal) G(eneral) A(array) C(omplex) H(ermitian) P(attern) S(ymmetric) I(nteger) K(kew) ***********************************************************************/ #define MM_MTX_STR "matrix" #define MM_ARRAY_STR "array" #define MM_DENSE_STR "array" #define MM_COORDINATE_STR "coordinate" #define MM_SPARSE_STR "coordinate" #define MM_COMPLEX_STR "complex" #define MM_REAL_STR "real" #define MM_INT_STR "integer" #define MM_GENERAL_STR "general" #define MM_SYMM_STR "symmetric" #define MM_HERM_STR "hermitian" #define MM_SKEW_STR "skew-symmetric" #define MM_PATTERN_STR "pattern" /* high level routines */ /* Matrix Market I/O library for ANSI C See https://math.nist.gov/MatrixMarket/ for details. */ #include <stdlib.h> #include <stdio.h> #include <string.h> #include <ctype.h> static char mm_buffer[MM_MAX_LINE_LENGTH]; int mm_read_unsymmetric_sparse(const char *fname, int *M_, int *N_, int *nz_, double **val_, int **I_, int **J_) { FILE *f; MM_typecode matcode; int M, N, nz; int i; double *val; int *ia, *ja; if ((f = fopen(fname, "r")) == NULL) return -1; if (mm_read_banner(f, &matcode) != 0) { printf("mm_read_unsymmetric: Could not process Matrix Market banner "); printf(" in file [%s]\n", fname); return -1; } if (!(mm_is_real(matcode) && mm_is_matrix(matcode) && mm_is_sparse(matcode))) { fprintf(stderr, "This application does not support "); fprintf(stderr, "Market Market type: [%s]\n", mm_typecode_to_str(matcode)); return -1; } /* find out size of sparse matrix: M, N, nz .... */ if (mm_read_mtx_crd_size(f, &M, &N, &nz) != 0) { fprintf(stderr, "read_unsymmetric_sparse(): could not parse matrix size.\n"); return -1; } *M_ = M; *N_ = N; *nz_ = nz; /* reseve memory for matrices */ ia = (int *)malloc(nz * sizeof(int)); ja = (int *)malloc(nz * sizeof(int)); val = (double *)malloc(nz * sizeof(double)); *val_ = val; *I_ = ia; *J_ = ja; /* NOTE: when reading in doubles, ANSI C requires the use of the "l" */ /* specifier as in "%lg", "%lf", "%le", otherwise errors will occur */ /* (ANSI C X3.159-1989, Sec. 4.9.6.2, p. 136 lines 13-15) */ for (i = 0; i < nz; i++) { if (fscanf(f, "%d %d %lg\n", &ia[i], &ja[i], &val[i]) != 3) { fprintf(stderr, "read_unsymmetric_sparse(): could not parse i, j and nonzero.\n"); return -1; } ia[i]--; /* adjust from 1-based to 0-based */ ja[i]--; } fclose(f); return 0; } int mm_is_valid(MM_typecode matcode) { if (!mm_is_matrix(matcode)) return 0; if (mm_is_dense(matcode) && mm_is_pattern(matcode)) return 0; if (mm_is_real(matcode) && mm_is_hermitian(matcode)) return 0; if (mm_is_pattern(matcode) && (mm_is_hermitian(matcode) || mm_is_skew(matcode))) return 0; return 1; } int mm_read_banner(FILE *f, MM_typecode *matcode) { char line[MM_MAX_LINE_LENGTH]; char banner[MM_MAX_TOKEN_LENGTH]; char mtx[MM_MAX_TOKEN_LENGTH]; char crd[MM_MAX_TOKEN_LENGTH]; char data_type[MM_MAX_TOKEN_LENGTH]; char storage_scheme[MM_MAX_TOKEN_LENGTH]; char *p; mm_clear_typecode(matcode); if (fgets(line, MM_MAX_LINE_LENGTH, f) == NULL) return MM_PREMATURE_EOF; if (sscanf(line, "%s %s %s %s %s", banner, mtx, crd, data_type, storage_scheme) != 5) return MM_PREMATURE_EOF; for (p = mtx; *p != '\0'; *p = tolower(*p), p++) ; /* convert to lower case */ for (p = crd; *p != '\0'; *p = tolower(*p), p++) ; for (p = data_type; *p != '\0'; *p = tolower(*p), p++) ; for (p = storage_scheme; *p != '\0'; *p = tolower(*p), p++) ; /* check for banner */ if (strncmp(banner, MatrixMarketBanner, strlen(MatrixMarketBanner)) != 0) return MM_NO_HEADER; /* first field should be "mtx" */ if (strcmp(mtx, MM_MTX_STR) != 0) return MM_UNSUPPORTED_TYPE; mm_set_matrix(matcode); /* second field describes whether this is a sparse matrix (in coordinate storgae) or a dense array */ if (strcmp(crd, MM_SPARSE_STR) == 0) mm_set_sparse(matcode); else if (strcmp(crd, MM_DENSE_STR) == 0) mm_set_dense(matcode); else return MM_UNSUPPORTED_TYPE; /* third field */ if (strcmp(data_type, MM_REAL_STR) == 0) mm_set_real(matcode); else if (strcmp(data_type, MM_COMPLEX_STR) == 0) mm_set_complex(matcode); else if (strcmp(data_type, MM_PATTERN_STR) == 0) mm_set_pattern(matcode); else if (strcmp(data_type, MM_INT_STR) == 0) mm_set_integer(matcode); else return MM_UNSUPPORTED_TYPE; /* fourth field */ if (strcmp(storage_scheme, MM_GENERAL_STR) == 0) mm_set_general(matcode); else if (strcmp(storage_scheme, MM_SYMM_STR) == 0) mm_set_symmetric(matcode); else if (strcmp(storage_scheme, MM_HERM_STR) == 0) mm_set_hermitian(matcode); else if (strcmp(storage_scheme, MM_SKEW_STR) == 0) mm_set_skew(matcode); else return MM_UNSUPPORTED_TYPE; return 0; } int mm_write_mtx_crd_size(FILE *f, int M, int N, int nz) { if (fprintf(f, "%d %d %d\n", M, N, nz) < 0) return MM_COULD_NOT_WRITE_FILE; else return 0; } int mm_read_mtx_crd_size(FILE *f, int *M, int *N, int *nz) { char line[MM_MAX_LINE_LENGTH]; int num_items_read; /* set return null parameter values, in case we exit with errors */ *M = *N = *nz = 0; /* now continue scanning until you reach the end-of-comments */ do { if (fgets(line, MM_MAX_LINE_LENGTH, f) == NULL) return MM_PREMATURE_EOF; } while (line[0] == '%'); /* line[] is either blank or has M,N, nz */ if (sscanf(line, "%d %d %d", M, N, nz) == 3) return 0; else do { num_items_read = fscanf(f, "%d %d %d", M, N, nz); if (num_items_read == EOF) return MM_PREMATURE_EOF; } while (num_items_read != 3); return 0; } int mm_read_mtx_array_size(FILE *f, int *M, int *N) { char line[MM_MAX_LINE_LENGTH]; int num_items_read; /* set return null parameter values, in case we exit with errors */ *M = *N = 0; /* now continue scanning until you reach the end-of-comments */ do { if (fgets(line, MM_MAX_LINE_LENGTH, f) == NULL) return MM_PREMATURE_EOF; } while (line[0] == '%'); /* line[] is either blank or has M,N, nz */ if (sscanf(line, "%d %d", M, N) == 2) return 0; else /* we have a blank line */ do { num_items_read = fscanf(f, "%d %d", M, N); if (num_items_read == EOF) return MM_PREMATURE_EOF; } while (num_items_read != 2); return 0; } int mm_write_mtx_array_size(FILE *f, int M, int N) { if (fprintf(f, "%d %d\n", M, N) < 0) return MM_COULD_NOT_WRITE_FILE; else return 0; } /*-------------------------------------------------------------------------*/ /******************************************************************/ /* use when ia[], ja[], and val[] are already allocated */ /******************************************************************/ int mm_read_mtx_crd_data(FILE *f, int M, int N, int nz, int ia[], int ja[], double val[], MM_typecode matcode) { int i; if (mm_is_complex(matcode)) { for (i = 0; i < nz; i++) if (fscanf(f, "%d %d %lg %lg", &ia[i], &ja[i], &val[2 * i], &val[2 * i + 1]) != 4) return MM_PREMATURE_EOF; } else if (mm_is_real(matcode)) { for (i = 0; i < nz; i++) { if (fscanf(f, "%d %d %lg\n", &ia[i], &ja[i], &val[i]) != 3) return MM_PREMATURE_EOF; } } else if (mm_is_pattern(matcode)) { for (i = 0; i < nz; i++) if (fscanf(f, "%d %d", &ia[i], &ja[i]) != 2) return MM_PREMATURE_EOF; } else return MM_UNSUPPORTED_TYPE; return 0; } int mm_read_mtx_crd_entry(FILE *f, int *ia, int *ja, double *real, double *imag, MM_typecode matcode) { if (mm_is_complex(matcode)) { if (fscanf(f, "%d %d %lg %lg", ia, ja, real, imag) != 4) return MM_PREMATURE_EOF; } else if (mm_is_real(matcode)) { if (fscanf(f, "%d %d %lg\n", ia, ja, real) != 3) return MM_PREMATURE_EOF; } else if (mm_is_pattern(matcode)) { if (fscanf(f, "%d %d", ia, ja) != 2) return MM_PREMATURE_EOF; } else return MM_UNSUPPORTED_TYPE; return 0; } /************************************************************************ mm_read_mtx_crd() fills M, N, nz, array of values, and return type code, e.g. 'MCRS' if matrix is complex, values[] is of size 2*nz, (nz pairs of real/imaginary values) ************************************************************************/ int mm_read_mtx_crd(char *fname, int *M, int *N, int *nz, int **ia, int **ja, double **val, MM_typecode *matcode) { int ret_code; FILE *f; if (strcmp(fname, "stdin") == 0) f = stdin; else if ((f = fopen(fname, "r")) == NULL) return MM_COULD_NOT_READ_FILE; if ((ret_code = mm_read_banner(f, matcode)) != 0) return ret_code; if (!(mm_is_valid(*matcode) && mm_is_sparse(*matcode) && mm_is_matrix(*matcode))) return MM_UNSUPPORTED_TYPE; if ((ret_code = mm_read_mtx_crd_size(f, M, N, nz)) != 0) return ret_code; *ia = (int *)malloc(*nz * sizeof(int)); *ja = (int *)malloc(*nz * sizeof(int)); *val = NULL; if (mm_is_complex(*matcode)) { *val = (double *)malloc(*nz * 2 * sizeof(double)); ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *ia, *ja, *val, *matcode); if (ret_code != 0) return ret_code; } else if (mm_is_real(*matcode)) { *val = (double *)malloc(*nz * sizeof(double)); ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *ia, *ja, *val, *matcode); if (ret_code != 0) return ret_code; } else if (mm_is_pattern(*matcode)) { ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *ia, *ja, *val, *matcode); if (ret_code != 0) return ret_code; } if (f != stdin) fclose(f); return 0; } int mm_write_banner(FILE *f, MM_typecode matcode) { char *str = mm_typecode_to_str(matcode); int ret_code; ret_code = fprintf(f, "%s %s\n", MatrixMarketBanner, str); if (ret_code < 0) return MM_COULD_NOT_WRITE_FILE; else return 0; } int mm_write_mtx_crd(char fname[], int M, int N, int nz, int ia[], int ja[], double val[], MM_typecode matcode) { FILE *f; int i; if (strcmp(fname, "stdout") == 0) f = stdout; else if ((f = fopen(fname, "w")) == NULL) return MM_COULD_NOT_WRITE_FILE; /* print banner followed by typecode */ fprintf(f, "%s ", MatrixMarketBanner); fprintf(f, "%s\n", mm_typecode_to_str(matcode)); /* print matrix sizes and nonzeros */ fprintf(f, "%d %d %d\n", M, N, nz); /* print values */ if (mm_is_pattern(matcode)) for (i = 0; i < nz; i++) fprintf(f, "%d %d\n", ia[i], ja[i]); else if (mm_is_real(matcode)) for (i = 0; i < nz; i++) fprintf(f, "%d %d %20.16g\n", ia[i], ja[i], val[i]); else if (mm_is_complex(matcode)) for (i = 0; i < nz; i++) fprintf(f, "%d %d %20.16g %20.16g\n", ia[i], ja[i], val[2 * i], val[2 * i + 1]); else { if (f != stdout) fclose(f); return MM_UNSUPPORTED_TYPE; } if (f != stdout) fclose(f); return 0; } char *mm_typecode_to_str(MM_typecode matcode) { const char *types[4]; /* check for MTX type */ if (mm_is_matrix(matcode)) types[0] = MM_MTX_STR; else return NULL; /* check for CRD or ARR matrix */ if (mm_is_sparse(matcode)) types[1] = MM_SPARSE_STR; else if (mm_is_dense(matcode)) types[1] = MM_DENSE_STR; else return NULL; /* check for element data type */ if (mm_is_real(matcode)) types[2] = MM_REAL_STR; else if (mm_is_complex(matcode)) types[2] = MM_COMPLEX_STR; else if (mm_is_pattern(matcode)) types[2] = MM_PATTERN_STR; else if (mm_is_integer(matcode)) types[2] = MM_INT_STR; else return NULL; /* check for symmetry type */ if (mm_is_general(matcode)) types[3] = MM_GENERAL_STR; else if (mm_is_symmetric(matcode)) types[3] = MM_SYMM_STR; else if (mm_is_hermitian(matcode)) types[3] = MM_HERM_STR; else if (mm_is_skew(matcode)) types[3] = MM_SKEW_STR; else return NULL; mm_buffer[0] = 0; snprintf(mm_buffer, sizeof(mm_buffer) / sizeof(mm_buffer[0]), "%s %s %s %s", types[0], types[1], types[2], types[3]); return mm_buffer; } #endif
#ifndef MM_LOADER_H
#define MM_LOADER_H
#include <petscmat.h>
#include "mmio.h"
PetscErrorCode MatCreateFromMTX(Mat *A, const char *filein, PetscBool aijonly)
{
MM_typecode matcode;
FILE *file;
PetscInt M, N, ninput;
PetscInt *ia, *ja;
PetscInt i, j, nz, *rownz;
PetscScalar *val;
PetscBool sametype, symmetric = PETSC_FALSE, skew = PETSC_FALSE;
/* Read in matrix */
PetscFunctionBeginUser;
PetscCall(PetscFOpen(PETSC_COMM_SELF, filein, "r", &file));
PetscCheck(mm_read_banner(file, &matcode) == 0, PETSC_COMM_SELF,
PETSC_ERR_ARG_WRONG, "Could not process Matrix Market banner.");
/* This is how one can screen matrix types if their application */
/* only supports a subset of the Matrix Market data types. */
PetscCheck(mm_is_matrix(matcode) && mm_is_sparse(matcode) &&
(mm_is_real(matcode) || mm_is_integer(matcode)), PETSC_COMM_WORLD,
PETSC_ERR_ARG_WRONG, "Input must be a sparse real or integer matrix. Market
Market type: [%s]", mm_typecode_to_str(matcode));
if (mm_is_symmetric(matcode))
symmetric = PETSC_TRUE;
if (mm_is_skew(matcode))
skew = PETSC_TRUE;
/* Find out size of sparse matrix .... */
PetscCheck(mm_read_mtx_crd_size(file, &M, &N, &nz) == 0, PETSC_COMM_SELF,
PETSC_ERR_ARG_WRONG, "Size of sparse matrix is wrong.");
/* Reserve memory for matrices */
PetscCall(PetscMalloc4(nz, &ia, nz, &ja, nz, &val, M, &rownz));
for (i = 0; i < M; i++)
rownz[i] = 0;
/* NOTE: when reading in doubles, ANSI C requires the use of the "l" */
/* specifier as in "%lg", "%lf", "%le", otherwise errors will occur */
/* (ANSI C X3.159-1989, Sec. 4.9.6.2, p. 136 lines 13-15) */
for (i = 0; i < nz; i++)
{
ninput = fscanf(file, "%d %d %lg\n", &ia[i], &ja[i], &val[i]);
PetscCheck(ninput >= 3, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED,
"Badly formatted input file");
if ((symmetric && aijonly) || skew)
{ /* transpose */
rownz[ia[i]]++;
rownz[ja[i]]++;
}
else
rownz[ia[i]]++;
}
PetscCall(PetscFClose(PETSC_COMM_SELF, file));
/* Create, preallocate, and then assemble the matrix */
PetscCall(MatCreate(PETSC_COMM_WORLD, A));
PetscCall(MatSetSizes(*A, PETSC_DECIDE, PETSC_DECIDE, M, N));
if (symmetric && !aijonly)
{
PetscCall(MatSetType(*A, MATSEQSBAIJ));
PetscCall(MatSetFromOptions(*A));
PetscCall(MatSeqSBAIJSetPreallocation(*A, 1, 0, rownz));
PetscCall(PetscObjectTypeCompare((PetscObject)(*A), MATSEQSBAIJ,
&sametype));
PetscCheck(sametype, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Only AIJ
and SBAIJ are supported. Your mattype is not supported");
}
else
{
PetscCall(MatSetType(*A, MATSEQAIJ));
PetscCall(MatSetFromOptions(*A));
PetscCall(MatSeqAIJSetPreallocation(*A, 0, rownz));
PetscCall(PetscObjectTypeCompare((PetscObject)(*A), MATSEQAIJ,
&sametype));
PetscCheck(sametype, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Only AIJ
and SBAIJ are supported. Your mattype is not supported");
}
/* Add values to the matrix, these correspond to lower triangular part for
symmetric or skew matrices */
for (j = 0; j < nz; j++)
PetscCall(MatSetValues(*A, 1, &ia[j], 1, &ja[j], &val[j],
INSERT_VALUES));
/* Add values to upper triangular part for some cases */
if (symmetric && aijonly)
{
/* MatrixMarket matrix stores symm matrix in lower triangular part.
Take its transpose */
for (j = 0; j < nz; j++)
PetscCall(MatSetValues(*A, 1, &ja[j], 1, &ia[j], &val[j],
INSERT_VALUES));
}
if (skew)
{
for (j = 0; j < nz; j++)
{
val[j] = -val[j];
PetscCall(MatSetValues(*A, 1, &ja[j], 1, &ia[j], &val[j],
INSERT_VALUES));
}
}
PetscCall(MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY));
PetscCall(MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY));
PetscCall(PetscFree4(ia, ja, val, rownz));
PetscFunctionReturn(PETSC_SUCCESS);
}
PetscErrorCode MatCreateFromMTX_1base(Mat *A, const char *filein, PetscBool
aijonly)
{
MM_typecode matcode;
FILE *file;
PetscInt M, N, ninput;
PetscInt *ia, *ja;
PetscInt i, j, nz, *rownz;
PetscScalar *val;
PetscBool sametype, symmetric = PETSC_FALSE, skew = PETSC_FALSE;
/* Read in matrix */
PetscFunctionBeginUser;
PetscCall(PetscFOpen(PETSC_COMM_WORLD, filein, "r", &file));
PetscCheck(mm_read_banner(file, &matcode) == 0, PETSC_COMM_WORLD,
PETSC_ERR_ARG_WRONG, "Could not process Matrix Market banner.");
/* This is how one can screen matrix types if their application */
/* only supports a subset of the Matrix Market data types. */
PetscCheck(mm_is_matrix(matcode) && mm_is_sparse(matcode) &&
(mm_is_real(matcode) || mm_is_integer(matcode)), PETSC_COMM_WORLD,
PETSC_ERR_ARG_WRONG, "Input must be a sparse real or integer matrix. Market
Market type: [%s]", mm_typecode_to_str(matcode));
if (mm_is_symmetric(matcode))
symmetric = PETSC_TRUE;
if (mm_is_skew(matcode))
skew = PETSC_TRUE;
/* Find out size of sparse matrix .... */
PetscCheck(mm_read_mtx_crd_size(file, &M, &N, &nz) == 0, PETSC_COMM_WORLD,
PETSC_ERR_ARG_WRONG, "Size of sparse matrix is wrong.");
/* Reserve memory for matrices */
PetscCall(PetscMalloc4(nz, &ia, nz, &ja, nz, &val, M, &rownz));
for (i = 0; i < M; i++)
rownz[i] = 0;
/* NOTE: when reading in doubles, ANSI C requires the use of the "l" */
/* specifier as in "%lg", "%lf", "%le", otherwise errors will occur */
/* (ANSI C X3.159-1989, Sec. 4.9.6.2, p. 136 lines 13-15) */
for (i = 0; i < nz; i++)
{
ninput = fscanf(file, "%d %d %lg\n", &ia[i], &ja[i], &val[i]);
PetscCheck(ninput >= 3, PETSC_COMM_WORLD, PETSC_ERR_FILE_UNEXPECTED,
"Badly formatted input file");
ia[i]--;
ja[i]--; /* adjust from 1-based to 0-based */
if ((symmetric && aijonly) || skew)
{ /* transpose */
rownz[ia[i]]++;
rownz[ja[i]]++;
}
else
rownz[ia[i]]++;
}
PetscCall(PetscFClose(PETSC_COMM_WORLD, file));
/* Create, preallocate, and then assemble the matrix */
PetscCall(MatCreate(PETSC_COMM_WORLD, A));
PetscCall(MatSetSizes(*A, PETSC_DECIDE, PETSC_DECIDE, M, N));
if (symmetric && !aijonly)
{
PetscCall(MatSetType(*A, MATSEQSBAIJ));
PetscCall(MatSetFromOptions(*A));
PetscCall(MatSeqSBAIJSetPreallocation(*A, 1, 0, rownz));
PetscCall(PetscObjectTypeCompare((PetscObject)(*A), MATSEQSBAIJ,
&sametype));
PetscCheck(sametype, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Only AIJ
and SBAIJ are supported. Your mattype is not supported");
}
else
{
PetscCall(MatSetType(*A, MATSEQAIJ));
PetscCall(MatSetFromOptions(*A));
PetscCall(MatSeqAIJSetPreallocation(*A, 0, rownz));
PetscCall(PetscObjectTypeCompare((PetscObject)(*A), MATSEQAIJ,
&sametype));
PetscCheck(sametype, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Only AIJ
and SBAIJ are supported. Your mattype is not supported");
}
/* Add values to the matrix, these correspond to lower triangular part for
symmetric or skew matrices */
for (j = 0; j < nz; j++)
PetscCall(MatSetValues(*A, 1, &ia[j], 1, &ja[j], &val[j],
INSERT_VALUES));
/* Add values to upper triangular part for some cases */
if (symmetric && aijonly)
{
/* MatrixMarket matrix stores symm matrix in lower triangular part.
Take its transpose */
for (j = 0; j < nz; j++)
PetscCall(MatSetValues(*A, 1, &ja[j], 1, &ia[j], &val[j],
INSERT_VALUES));
}
if (skew)
{
for (j = 0; j < nz; j++)
{
val[j] = -val[j];
PetscCall(MatSetValues(*A, 1, &ja[j], 1, &ia[j], &val[j],
INSERT_VALUES));
}
}
PetscCall(MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY));
PetscCall(MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY));
PetscCall(PetscFree4(ia, ja, val, rownz));
PetscFunctionReturn(PETSC_SUCCESS);
}
PetscErrorCode VecCreateFromRHS(Vec *A, const char *filein)
{
FILE *file;
PetscInt nz;
PetscScalar *val;
PetscInt *idx;
PetscInt i;
PetscCall(PetscFOpen(PETSC_COMM_WORLD, filein, "r", &file));
fscanf(file, "%d", &nz);
PetscCall(PetscMalloc1(nz, &val));
PetscCall(PetscMalloc1(nz, &idx));
for (i = 0; i < nz; i++)
{
fscanf(file, "%lg", &val[i]);
idx[i] = i;
}
printf("read done\n");
PetscCall(VecCreate(PETSC_COMM_WORLD, A));
printf("VecSetSizes\n");
PetscCall(VecSetSizes(*A, PETSC_DECIDE, nz));
PetscCall(VecSetType(*A, VECMPI));
// printf("VecSet\n");
// PetscCall(VecSet(*A, 1.0));
printf("VecSetValues\n");
PetscCall(VecSetValues(*A, nz, idx, val, INSERT_VALUES));
printf("VecAssemblyBegin\n");
PetscCall(VecAssemblyBegin(*A));
printf("VecAssemblyEnd\n");
PetscCall(VecAssemblyEnd(*A));
printf("PetscFree\n");
PetscCall(PetscFree(val));
PetscFunctionReturn(PETSC_SUCCESS);
}
PetscErrorCode load_binary_mat(Mat *A, char *file)
{
PetscErrorCode ierr;
PetscViewer viewer;
// ierr = PetscPrintf(PETSC_COMM_WORLD, "Read Matrix in binary
format...\n");
// CHKERRQ(ierr);
ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ,
&viewer);
CHKERRQ(ierr);
ierr = MatCreate(PETSC_COMM_WORLD, A);
CHKERRQ(ierr);
ierr = MatSetType(*A, MATAIJ);
CHKERRQ(ierr);
ierr = MatLoad(*A, viewer);
CHKERRQ(ierr);
ierr = PetscViewerDestroy(&viewer);
CHKERRQ(ierr);
return ierr;
}
PetscErrorCode load_binary_rhs(Vec *b, char *file)
{
PetscErrorCode ierr;
PetscViewer viewer;
// ierr = PetscPrintf(PETSC_COMM_WORLD, "Read rhs in binary format...\n");
// CHKERRQ(ierr);
ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ,
&viewer);
CHKERRQ(ierr);
ierr = VecCreate(PETSC_COMM_WORLD, b);
CHKERRQ(ierr);
ierr = VecSetType(*b, VECMPI);
CHKERRQ(ierr);
ierr = VecLoad(*b, viewer);
CHKERRQ(ierr);
ierr = PetscViewerDestroy(&viewer);
CHKERRQ(ierr);
return ierr;
}
// double vec2norm(double *x, int n)
// {
// double sum = 0.0;
// for (int i = 0; i < n; i++)
// sum += x[i] * x[i];
// return sqrt(sum);
// }
// calculate_type check_calculate_x(calculate_type **x_ref, calculate_type
**cal, int N, int nrhs)
// {
// calculate_type e = 0.0;
// for (int r = 0; r < nrhs; r++)
// {
// calculate_type *ref_answer = x_ref[r];
// calculate_type *answer = cal[r];
// for (int i = 0; i < N; i++)
// {
// answer[i] -= ref_answer[i];
// }
// for (int i = 0; i < N; i++)
// {
// if (answer[i] > ERROR || answer[i] < -ERROR)
// { // printf("ROW %d ANSWER %lf ERROR %lf\n", i, ref_answer[i],
answer[i]);
// }
// }
// e += vec2norm(answer, N) / vec2norm(ref_answer, N);
// }
// return e;
// }
#endif
slurm-4803840.out
Description: Binary data
task.sh
Description: Bourne shell script
