Commit 7f93d0b7 authored by tempy@maggie's avatar tempy@maggie
Browse files

Common functions moved to one place

parent b6e5124a
#ifndef CSR_H
#define CSR_H
typedef int csi;
typedef double csv;
csv zero = 0.0;
#define CS_MAX(a,b) (((a) > (b)) ? (a) : (b))
typedef struct csr_t {
csi m;
csi n;
csi nzmax;
csi nr;
csi* r;
csi* p;
csi* j;
csv* x;
} csr;
/* method declarations {{{*/
csr *csr_spfree(csr *A);
/*}}}*/
/*csr util{{{*/
/* free workspace and return a sparse matrix result */
csr *csr_done(csr *C, void *w, void *x, csi ok) {
return(ok ? C : csr_spfree(C)); /* return result if OK, else free it */
}
/* wrapper for free */
void *cs_free(void *p) {
if (p)
free(p); /* free p if it is not already NULL */
return(NULL); /* return NULL to simplify the use of cs_free */
}
/* wrapper for realloc */
void *csr_realloc(void *p, csi n, size_t size, csi *ok) {
void *pnew = NULL;
pnew = realloc(p, CS_MAX(n, 1) * size); /* realloc the block */
*ok = (pnew != NULL); /* realloc fails if pnew is NULL */
if (pnew == NULL) {
printf("%d:reallocation failed, pnew is NULL\n", __LINE__);
}
return((*ok) ? pnew : p); /* return original p if failure */
}
/* wrapper for realloc */
void *cs_realloc(void *p, csi n, size_t size, csi *ok) {
void *pnew = NULL;
pnew = realloc(p, CS_MAX(n, 1) * size); /* realloc the block */
*ok = (pnew != NULL); /* realloc fails if pnew is NULL */
if (pnew == NULL) {
printf("reallocation failed\n");
}
return((*ok) ? pnew : p); /* return original p if failure */
}
/* change the max # of entries sparse matrix */
csi csr_sprealloc(csr *A, csi nzmax) {
csi ok, oki = 0, okj = 1, okx = 1;
if (!A)
return(0);
if (nzmax <= 0)
nzmax = A->p[A->m];
A->j = (int*)csr_realloc(A->j, nzmax, sizeof(csi), &oki);
if (A->x)
A->x = (csv*)csr_realloc(A->x, nzmax, sizeof(csv), &okx);
ok = (oki && okj && okx);
if (ok)
A->nzmax = nzmax;
return(ok);
}
/* free a sparse matrix */
csr *csr_spfree(csr *A) {
if (!A)
return(NULL); /* do nothing if A already NULL */
cs_free(A->p);
A->p = NULL;
cs_free(A->j);
A->j = NULL;
cs_free(A->x);
A->x = NULL;
cs_free(A->r);
A->r = NULL;
cs_free(A); /* free the cs struct and return NULL */
return NULL;
}
/* allocate a sparse matrix (triplet form or compressed-ROW form) */
csr *csr_spalloc(csi m, csi n, csi nzmax, int values, int triplet, csv f) {
csr* A = (csr*)calloc(1, sizeof(csr)); /* allocate the cs struct */
if (!A) {
perror("sparse allocation failed");
return(NULL); /* out of memory */
}
A->m = m; /* define dimensions and nzmax */
A->n = n;
A->nzmax = nzmax = CS_MAX(nzmax, 0);
A->nr = 0; // number of nonzero rows
A->p = (csi*)calloc(m + 2, sizeof(csi));
A->j = (csi*)calloc(CS_MAX(nzmax,1), sizeof(csi));
A->x = (csv*)calloc(CS_MAX(nzmax,1), sizeof(csv));
return((!A->p || !A->j || !A->x) ? csr_spfree(A) : A);
}/*}}}*/
/** Multiply two sparse matrices which are stored in CSR format. MKL is used */
csr *csr_multiply(csi Am, csi An, csi Anzmax, const csi* Ap, const csi* Aj, const csv* Ax, csi Bm, csi Bn, csi Bnzmax, const csi* Bp, const csi* Bj, const csv* Bx, long* nummult, csi* xb, csv* x) { /*{{{*/
csv tf = 0;
csi p, jp, j, kp, k, i, nz = 0, anz, *Cp, *Cj, m, n,
bnz, values = 1;
csv *Cx;
csr *C;
if (An != Bm)
return(NULL);
if (Anzmax == 0 || Bnzmax == 0) {
C = csr_spalloc(Am, Bn, 0, values, 0, tf);
return C;
}
m = Am;
anz = Ap[Am];
n = Bn;
bnz = Bp[Bm];
for(i = 0; i < n; i++) xb[i] = 0;
for(i = 0; i < n; i++)
xb[i] = 0;
values = (Ax != NULL) && (Bx != NULL);
csi tnz = (anz + bnz) * 2;
C = csr_spalloc(m, n, tnz, values, 0, tf); /* allocate result */
if (!C || !xb || (values && !x))
return (csr_done(C, xb, x, 0));
Cp = C->p;
for (i = 0; i < m; i++) {
if ( ( (nz + n) > C->nzmax ) ) {
if(!csr_sprealloc(C, (2 * (C->nzmax) + n) ) ) {
return (csr_done(C, xb, x, 0)); // out of memory
} else {
}
}
Cj = C->j;
Cx = C->x; /* C->j and C->x may be reallocated */
Cp[i] = nz; /* row i of C starts here */
for (jp = Ap[i]; jp < Ap[i + 1]; jp++) {
j = Aj[jp];
for (kp = Bp[j]; kp < Bp[j + 1]; kp++) {
k = Bj[kp]; /* B(i,j) is nonzero */
if (xb[k] != i + 1) {
xb[k] = i + 1; /* i is new entry in column j */
Cj[nz++] = k; /* add i to pattern of C(:,j) */
if (x) {
x[k] = Ax[jp] * Bx[kp]; /* x(i) = beta*A(i,j) */
(*nummult)++;
}
} else if (x) {
x[k] += (Ax[jp] * Bx[kp]); /* i exists in C(:,j) already */
(*nummult)++;
}
}
}
if (values)
for (p = Cp[i]; p < nz; p++)
Cx[p] = x[Cj[p]];
}
Cp[m] = nz; /* finalize the last row of C */
csr_sprealloc(C, 0); /* remove extra space from C */
xb = NULL;
x = NULL;
return C;
}/*}}}*/
#endif
......@@ -6,14 +6,133 @@
*
*/
#ifndef MM_IO_H
#define MM_IO_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 fucntions ***************************/
#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 fucntions ***************************/
#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
ojbect 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 */
int mm_write_mtx_crd(char fname[], int M, int N, int nz, int I[], int J[],
double val[], MM_typecode matcode);
int mm_read_mtx_crd_data(FILE *f, int M, int N, int nz, int I[], int J[],
double val[], MM_typecode matcode);
int mm_read_mtx_crd_entry(FILE *f, int *I, int *J, double *real, double *img,
MM_typecode matcode);
int mm_read_unsymmetric_sparse(const char *fname, int *M_, int *N_, int *nz_,
double **val_, int **I_, int **J_);
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <ctype.h>
#include "mmio.h"
int mm_read_unsymmetric_sparse(const char *fname, int *M_, int *N_, int *nz_,
double **val_, int **I_, int **J_)
{
......@@ -509,3 +628,66 @@ char *mm_typecode_to_str(MM_typecode matcode)
return mm_strdup(buffer);
}
/** Read Matrix Market file into COO matrix */
void read_mm(char* strpath, int* pM, int* pN, int* prealnnz, int** pI, int** pJ, double** pval){ /*{{{*/
/*
* taken from Matrix Market I/O library for ANSI C
*
* See http://math.nist.gov/MatrixMarket for details.
*
*
*/
int i, M, N, nz, *I, *J;
double* val;
int ret_code;
MM_typecode matcode;
FILE* f;
if ((f = fopen(strpath, "r")) == NULL) {fprintf(stderr, "Input matrix file %s cannot be opened to read.", strpath);exit(1);}
/* READ MATRIX */
if (mm_read_banner(f, &matcode) != 0) {
printf("Could not process Matrix Market banner.\n"); exit(1);
}
/* This is how one can screen matrix types if their application */
/* only supports a subset of the Matrix Market data types. */
if (mm_is_complex(matcode) && mm_is_matrix(matcode) && mm_is_sparse(matcode) ) {
printf("Sorry, this application does not support ");
printf("Market Market type: [%s]\n", mm_typecode_to_str(matcode));exit(1);
}
/* find out size of sparse matrix .... */
if ((ret_code = mm_read_mtx_crd_size(f, &M, &N, &nz)) !=0) exit(1);
/* reseve memory for matrices */
I = (int *) malloc((mm_is_symmetric(matcode) ? 2*nz : nz) * sizeof(int));
J = (int *) malloc((mm_is_symmetric(matcode) ? 2*nz : nz) * sizeof(int));
val = (double *) malloc((mm_is_symmetric(matcode) ? 2*nz : nz) * sizeof(double));
*pI = I;
*pJ = J;
*pval = val;
/* 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) */
int realnnz = 0;
for (i=0; i<nz; i++) {
if(mm_is_pattern(matcode)) {
fscanf(f, "%d %d\n", &I[realnnz], &J[realnnz]);
val[realnnz] = 1.0;
}
else
fscanf(f, "%d %d %lg\n", &I[realnnz], &J[realnnz], &val[realnnz]);
I[realnnz]--; /* adjust from 1-based to 0-based */
J[realnnz]--;
if(mm_is_symmetric(matcode) && I[realnnz] != J[realnnz]) {
I[realnnz+1] = J[realnnz];
J[realnnz+1] = I[realnnz];
val[realnnz+1] = val[realnnz];
realnnz++;
}
realnnz++;
}
if (f !=stdin) fclose(f);
*pM = M;
*pN = N;
*prealnnz = realnnz;
} /* ENDOF read_mm }}}*/
#endif
#ifndef UTIL_H
#define UTIL_H
#define OPTION_NOPRINT_MATRICES 0
#define OPTION_PRINT_MATRICES 1
MIC_TARGET int option_print_matrices = 0;
/** Converts COO matrix to CSR matrix */
void coo_to_csr(int m, int nnz, int* I, int* J, double* val, MKL_INT* AI, MKL_INT* AJ, double* Aval) { /*{{{*/
MKL_INT info = 0;
MKL_INT job[8];
job[1]=1; // one based indexing in csr
job[2]=0; // zero based indexing in coo
job[3]=2; // I don't know
job[4]=nnz; // nnz
job[0]=1; // coo to csr
job[5]=0; // Acsr and AJR allocated by user
mkl_dcsrcoo (job,&m, Aval, AJ, AI, &nnz, val, I, J, &info);
} /* ENDOF coo_to_csr }}}*/
/** Prints matrix in CSR format */
void MIC_TARGET printmm_one(int m, double* Aval, int* AJ, int* AI){ //{{{
if (option_print_matrices == OPTION_NOPRINT_MATRICES)
return;
int i;
for(i = 0; i < m; i++) {
printf("%d: ", i+1);
int j;
for(j = AI[i]-1; j < AI[i+1]-1; j++) {
printf("%d:%g ", AJ[j], Aval[j]);
}
printf("\n");
}
printf("\n");
}//}}}
/** Writes matrix in CSR format in to a file using Matrix Market format */
void MIC_TARGET printfilemm_one(char* file, int m, int n, double* Aval, int* AJ, int* AI){//{{{
FILE* f = fopen(file, "w");
if(f == NULL){
printf("%s %s %d: %s cannot be opened to write matrix\n", __FILE__, __PRETTY_FUNCTION__, __LINE__, file);
exit(1);
}
int i;
fprintf(f, "%%%%MatrixMarket matrix coordinate real general\n");
fprintf(f, "%d %d %d\n", m, n, AI[m]-1);
for(i = 0; i < m; i++) {
int j;
for(j = AI[i]-1; j < AI[i+1]-1; j++) {
fprintf(f, "%d %d %g\n", i+1, AJ[j], Aval[j]);
}
}
fclose(f);
}//}}}
#endif
......@@ -47,7 +47,7 @@ set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${CXX_FLAGS}")
if (MKL_FOUND)
include_directories(${MKL_INCLUDE_DIR})
link_directories(${MKL_LIBRARY_DIR})
add_executable(${NAME} mklspgemm.c mmio.c)
add_executable(${NAME} mklspgemm.c)
target_link_libraries(${NAME} mkl_intel_lp64 mkl_sequential mkl_core pthread m)
install(TARGETS ${NAME} DESTINATION bin)
else ()
......
......@@ -3,424 +3,132 @@
#include <assert.h>
#include <math.h>
#include <mkl.h>
#include "mmio.h"
#define MIC_TARGET
#include "../../common/mmio.h"
#include "../../common/util.h"
#include "../../common/csr.h"
/*}}}*/
/* global varibles and defines {{{*/
int nrepeat = 100;
int option_print_matrices = 0;
#define OPTION_NOPRINT_MATRICES 0
#define OPTION_PRINT_MATRICES 1
char transa = 'n';
typedef int csi;
typedef double csv;
csv zero = 0.0;
#define CS_MAX(a,b) (((a) > (b)) ? (a) : (b))
typedef struct csr_t {
csi m;
csi n;
csi nzmax;
csi nr;
csi* r;
csi* p;
csi* j;
csv* x;
} csr;
/*}}}*/
/* method declarations {{{*/
int main(int argc, char* argv[]);
void printmm_one(int m, double* Aval, int* AJ, int* AI);
void printfilemm_one(char* file, int m, int n, double* Aval, int* AJ, int* AI);
void printmm_zero(int m, double* Aval, int* AJ, int* AI);
csr *csr_spfree(csr *A);
/*}}}*/
/*csr util{{{*/
/* free workspace and return a sparse matrix result */
csr *csr_done(csr *C, void *w, void *x, csi ok) {
return(ok ? C : csr_spfree(C)); /* return result if OK, else free it */
}
/* wrapper for free */
void *cs_free(void *p) {
if (p)
free(p); /* free p if it is not already NULL */
return(NULL); /* return NULL to simplify the use of cs_free */
}
/* wrapper for realloc */
void *csr_realloc(void *p, csi n, size_t size, csi *ok) {
void *pnew = NULL;
pnew = realloc(p, CS_MAX(n, 1) * size); /* realloc the block */
*ok = (pnew != NULL); /* realloc fails if pnew is NULL */
if (pnew == NULL) {
printf("%d:reallocation failed, pnew is NULL\n", __LINE__);
}
//printf("%s:%d: n=%d ok=%d\n", __FUNCTION__, __LINE__, n, *ok);
return((*ok) ? pnew : p); /* return original p if failure */
}
/* wrapper for realloc */
void *cs_realloc(void *p, csi n, size_t size, csi *ok) {
void *pnew = NULL;
pnew = realloc(p, CS_MAX(n, 1) * size); /* realloc the block */
*ok = (pnew != NULL); /* realloc fails if pnew is NULL */
if (pnew == NULL) {
printf("reallocation failed\n");
}
return((*ok) ? pnew : p); /* return original p if failure */
}
/* change the max # of entries sparse matrix */
csi csr_sprealloc(csr *A, csi nzmax) {
csi ok, oki = 0, okj = 1, okx = 1;
if (!A)
return(0);
if (nzmax <= 0)
nzmax = A->p[A->m];
A->j = (int*)csr_realloc(A->j, nzmax, sizeof(csi), &oki);
if (A->x)
A->x = (csv*)csr_realloc(A->x, nzmax, sizeof(csv), &okx);
ok = (oki && okj && okx);
if (ok)
A->nzmax = nzmax;
return(ok);
}
/* free a sparse matrix */
csr *csr_spfree(csr *A) {
if (!A)
return(NULL); /* do nothing if A already NULL */
cs_free(A->p);
A->p = NULL;
cs_free(A->j);
A->j = NULL;
cs_free(A->x);
A->x = NULL;
cs_free(A->r);
A->r = NULL;
cs_free(A); /* free the cs struct and return NULL */
return NULL;
}
/* allocate a sparse matrix (triplet form or compressed-ROW form) */
csr *csr_spalloc(csi m, csi n, csi nzmax, int values, int triplet, csv f) {
csr* A = (csr*)calloc(1, sizeof(csr)); /* allocate the cs struct */
if (!A) {
perror("sparse allocation failed");
return(NULL); /* out of memory */
}
A->m = m; /* define dimensions and nzmax */
A->n = n;
A->nzmax = nzmax = CS_MAX(nzmax, 0);
A->nr = 0; // number of nonzero rows
A->p = (csi*)calloc(m + 2, sizeof(csi));
A->j = (csi*)calloc(CS_MAX(nzmax,1), sizeof(csi));
A->x = (csv*)calloc(CS_MAX(nzmax,1), sizeof(csv));
return((!A->p || !A->j || !A->x) ? csr_spfree(A) : A);
}/*}}}*/
/** Multiply two sparse matrices which are stored in CSR format. MKL is used */
csr *csr_multiply(csi Am, csi An, csi Anzmax, const csi* Ap, const csi* Aj, const csv* Ax, csi Bm, csi Bn, csi Bnzmax, const csi* Bp, const csi* Bj, const csv* Bx, long* nummult, csi* xb, csv* x) { /*{{{*/
csv tf = 0;
csi p, jp, j, kp, k, i, nz = 0, anz, *Cp, *Cj, m, n,
bnz, values = 1;
csv *Cx;
csr *C;
if (An != Bm)
return(NULL);
if (Anzmax == 0 || Bnzmax == 0) {
C = csr_spalloc(Am, Bn, 0, values, 0, tf);
return C;
}
m = Am;
anz = Ap[Am];
n = Bn;
bnz = Bp[Bm];
for(i = 0; i < n; i++) xb[i] = 0;
for(i = 0; i < n; i++)
xb[i] = 0;
values = (Ax != NULL) && (Bx != NULL);
csi tnz = (anz + bnz) * 2;
C = csr_spalloc(m, n, tnz, values, 0, tf); /* allocate result */
if (!C || !xb || (values && !x))