Skip to content
GitLab
Commit 3796eb48 authored by kadir-hpi7's avatar kadir-hpi7
Browse files

I have completed mkl_shmem

parent 51e0a644
Expand all Hide whitespace changes
Inline Side-by-side
spgemm/mkl_shmem/COMMAND.HISTORY.txt deleted 100644 → 0
View file @ 51e0a644
make;./mklspmm ../mtx/cp2k-h2o-e6.mtx 30 100 null
spgemm/mkl_shmem/Makefile
View file @ 3796eb48
CFLAGS:=-O0 -g3 #FIXME
#CFLAGS:=-O1
#CFLAGS:=-O0 -g3
CFLAGS:=-O3
all: mklspmm
all: mklspgemm
#all: mklspmv mklspmm
mklspmv: mklspmv.o mmio.o
icc mklspmv.o mmio.o -o $@ -mkl -lm -offload-option,mic,ld,"-rpath /opt/intel/mic/myo/lib/ -rpath /opt/intel/mic/coi/device-linux-release/lib/"
mklspmm: mklspmm.o mmio.o
icc mklspmm.o mmio.o -o $@ -mkl -lm
mklspgemm: mklspgemm.o mmio.o
icc mklspgemm.o mmio.o -o $@ -mkl -lm
.c.o:
icc $(CFLAGS) $(OPTION) -o $@ -c $<
run:
./mv ../mtx/tinyA.mtx
clean:
@rm mmio.o mklspmv mklspmm mklspmm.o mklspmv.o
@rm mmio.o mklspgemm mklspgemm.o
spgemm/mkl_shmem/QUICKSTART.txt
View file @ 3796eb48
make;./mklspmm ../mtx/smallA.mtx ../mtx/smallA.mtx 30 100 -1 null null perm5.txt PRINTMATRICES
make
./mklspgemm test.mtx test.mtx out.mtx 2 PRINT_YES
spgemm/mkl_shmem/common_func.c deleted 100644 → 0
View file @ 51e0a644
This diff is collapsed.
spgemm/mkl_shmem/commonmm.h deleted 100644 → 0
View file @ 51e0a644
#ifndef __COMMONMM_H__
#define __COMMONMM_H_
#define FIELD_LENGTH 128
//char filename[FIELD_LENGTH];
enum spmv_target { use_cpu, use_mkl, use_mic, use_mkl_mic, use_mic_rowcyclic };
char *target_str[] = { "CPU", "MKL", "MIC", "MKL_MIC", "MIC_ROWCYCLIC" };
enum enum_datastructure { DS_ICSR, DS_CSR };
char *str_datastructure[] = { "ICSR", "CSR" };
char *str_algorithm[] = { "OUTER", "INNER" };
__declspec(target(mic)) int datastructure;
typedef int csi;
typedef double csv;
__declspec(target(mic)) csv zero = 0.0;
__declspec(target(mic)) int use_sparse_c;
__declspec(target(mic)) int size_dacc;
__declspec(target(mic)) csv* dacc;
typedef struct csr_t {
csi m;
csi n;
csi nzmax;
csi nr;
csi* r;
csi* p;
csi* j;
csv* x;
int numInMats; // number of input matrices stored in this struct
int numRows;
int nnonempty_rows;
int numCols;
int totalnnz;
// int totalnnzmax;
int* nItems; // nnz of each matrix
int* pirow;
int size_pirow;
int* pnr;
int* h_rowDelimiters;
int size_h_rowDelimiters;
int* h_cols;
csv* h_val;
char filename[FIELD_LENGTH];
} csr;
#endif
spgemm/mkl_shmem/example_read.c deleted 100644 → 0
View file @ 51e0a644
/*
* Matrix Market I/O example program
*
* Read a real (non-complex) sparse matrix from a Matrix Market (v. 2.0) file.
* and copies it to stdout. This porgram does nothing useful, but
* illustrates common usage of the Matrix Matrix I/O routines.
* (See http://math.nist.gov/MatrixMarket for details.)
*
* Usage: a.out [filename] > output
*
*
* NOTES:
*
* 1) Matrix Market files are always 1-based, i.e. the index of the first
* element of a matrix is (1,1), not (0,0) as in C. ADJUST THESE
* OFFSETS ACCORDINGLY offsets accordingly when reading and writing
* to files.
*
* 2) ANSI C requires one to use the "l" format modifier when reading
* double precision floating point numbers in scanf() and
* its variants. For example, use "%lf", "%lg", or "%le"
* when reading doubles, otherwise errors will occur.
*/
#include <stdio.h>
#include <stdlib.h>
#include "mmio.h"
int main(int argc, char *argv[])
{
int ret_code;
MM_typecode matcode;
FILE *f;
int M, N, nz;
int i, *I, *J;
double *val;
if (argc < 2)
{
fprintf(stderr, "Usage: %s [martix-market-filename]\n", argv[0]);
exit(1);
}
else
{
if ((f = fopen(argv[1], "r")) == NULL)
exit(1);
}
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(nz * sizeof(int));
J = (int *) malloc(nz * sizeof(int));
val = (double *) malloc(nz * sizeof(double));
/* 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++)
{
fscanf(f, "%d %d %lg\n", &I[i], &J[i], &val[i]);
I[i]--; /* adjust from 1-based to 0-based */
J[i]--;
}
if (f !=stdin) fclose(f);
/************************/
/* now write out matrix */
/************************/
mm_write_banner(stdout, matcode);
mm_write_mtx_crd_size(stdout, M, N, nz);
for (i=0; i<nz; i++)
fprintf(stdout, "%d %d %20.19g\n", I[i]+1, J[i]+1, val[i]);
return 0;
}
spgemm/mkl_shmem/example_write.c deleted 100644 → 0
View file @ 51e0a644
/*
* Matrix Market I/O example program
*
* Create a small sparse, coordinate matrix and print it out
* in Matrix Market (v. 2.0) format to standard output.
*
* (See http://math.nist.gov/MatrixMarket for details.)
*
*/
#include <stdio.h>
#include <stdlib.h>
#include "mmio.h"
#define nz 4
#define M 10
#define N 10
int main()
{
MM_typecode matcode;
int I[nz] = { 0, 4, 2, 8 };
int J[nz] = { 3, 8, 7, 5 };
double val[nz] = {1.1, 2.2, 3.2, 4.4};
int i;
mm_initialize_typecode(&matcode);
mm_set_matrix(&matcode);
mm_set_coordinate(&matcode);
mm_set_real(&matcode);
mm_write_banner(stdout, matcode);
mm_write_mtx_crd_size(stdout, M, N, nz);
/* NOTE: matrix market files use 1-based indices, i.e. first element
of a vector has index 1, not 0. */
for (i=0; i<nz; i++)
fprintf(stdout, "%d %d %10.3g\n", I[i]+1, J[i]+1, val[i]);
return 0;
}
spgemm/mkl_shmem/mkl_example.h deleted 100644 → 0
View file @ 51e0a644
/*******************************************************************************
! Copyright(C) 1999-2013 Intel Corporation. All Rights Reserved.
!
! The source code, information and material ("Material") contained herein is
! owned by Intel Corporation or its suppliers or licensors, and title to such
! Material remains with Intel Corporation or its suppliers or licensors. The
! Material contains proprietary information of Intel or its suppliers and
! licensors. The Material is protected by worldwide copyright laws and treaty
! provisions. No part of the Material may be used, copied, reproduced,
! modified, published, uploaded, posted, transmitted, distributed or disclosed
! in any way without Intel's prior express written permission. No license
! under any patent, copyright or other intellectual property rights in the
! Material is granted to or conferred upon you, either expressly, by
! implication, inducement, estoppel or otherwise. Any license under such
! intellectual property rights must be express and approved by Intel in
! writing.
!
! *Third Party trademarks are the property of their respective owners.
!
! Unless otherwise agreed by Intel in writing, you may not remove or alter
! this notice or any other notice embedded in Materials by Intel or Intel's
! suppliers or licensors in any way.
!
!*******************************************************************************
! Content:
! Intel(R) Math Kernel Library (MKL) examples interface
!******************************************************************************/
#ifndef _MKL_EXAMPLE_H_
#define _MKL_EXAMPLE_H_
#include "mkl_types.h"
#include "mkl_cblas.h"
#include <stdio.h>
#ifdef __cplusplus
extern "C" {
#endif /* __cplusplus */
#define COMMENTS ':'
#define MAX_STRING_LEN 512
#define MIN(a,b) (((a) > (b)) ? (b) : (a))
#define MAX(a,b) (((a) > (b)) ? (a) : (b))
#define GENERAL_MATRIX 0
#define UPPER_MATRIX 1
#define LOWER_MATRIX -1
#define FULLPRINT 0
#define SHORTPRINT 1
#ifdef MKL_ILP64
#define INT_FORMAT "%lld"
#else
#define INT_FORMAT "%d"
#endif
/*
* ===========================================
* Prototypes for example program functions
* ===========================================
*/
int GetIntegerParameters(FILE*, ...);
int GetScalarsD(FILE*, ...);
int GetScalarsS(FILE*, ...);
int GetScalarsC(FILE*, ...);
int GetScalarsZ(FILE*, ...);
int GetCblasCharParameters(FILE *in_file, ...);
MKL_INT MaxValue(MKL_INT, MKL_INT*);
MKL_INT GetVectorI(FILE*, MKL_INT, MKL_INT*);
MKL_INT GetVectorS(FILE*, MKL_INT, float*, MKL_INT);
MKL_INT GetVectorC(FILE*, MKL_INT, MKL_Complex8*, MKL_INT);
MKL_INT GetVectorD(FILE*, MKL_INT, double*, MKL_INT);
MKL_INT GetVectorZ(FILE*, MKL_INT, MKL_Complex16*, MKL_INT);
MKL_INT GetArrayS(FILE*, CBLAS_ORDER*, MKL_INT, MKL_INT*, MKL_INT*, float*, MKL_INT*);
MKL_INT GetArrayD(FILE*, CBLAS_ORDER*, MKL_INT, MKL_INT*, MKL_INT*, double*, MKL_INT*);
MKL_INT GetArrayC(FILE*, CBLAS_ORDER*, MKL_INT, MKL_INT*, MKL_INT*, MKL_Complex8*, MKL_INT*);
MKL_INT GetArrayZ(FILE*, CBLAS_ORDER*, MKL_INT, MKL_INT*, MKL_INT*, MKL_Complex16*, MKL_INT*);
MKL_INT GetBandArrayS(FILE*, CBLAS_ORDER*, MKL_INT, MKL_INT, MKL_INT, MKL_INT, float*, MKL_INT);
MKL_INT GetBandArrayD(FILE*, CBLAS_ORDER*, MKL_INT, MKL_INT, MKL_INT, MKL_INT, double*, MKL_INT);
MKL_INT GetBandArrayC(FILE*, CBLAS_ORDER*, MKL_INT, MKL_INT, MKL_INT, MKL_INT, MKL_Complex8*, MKL_INT);
MKL_INT GetBandArrayZ(FILE*, CBLAS_ORDER*, MKL_INT, MKL_INT, MKL_INT, MKL_INT, MKL_Complex16*, MKL_INT);
MKL_INT GetValuesI(FILE *, MKL_INT *, MKL_INT, MKL_INT);
MKL_INT GetValuesC(FILE *, MKL_Complex8*, MKL_INT, MKL_INT, MKL_INT);
MKL_INT GetValuesZ(FILE *, MKL_Complex16*, MKL_INT, MKL_INT, MKL_INT);
MKL_INT GetValuesD(FILE *, double*, MKL_INT, MKL_INT, MKL_INT);
MKL_INT GetValuesS(FILE *, float*, MKL_INT, MKL_INT, MKL_INT);
void PrintParameters( char*, ... );
void PrintTrans(MKL_INT, CBLAS_TRANSPOSE*, CBLAS_TRANSPOSE*);
void PrintOrder(CBLAS_ORDER*);
void PrintVectorI(MKL_INT, MKL_INT*, char*);
void PrintVectorS(int, MKL_INT, float*, MKL_INT, char*);
void PrintVectorC(int, MKL_INT, MKL_Complex8*, MKL_INT, char*);
void PrintVectorD(int, MKL_INT, double*, MKL_INT, char*);
void PrintVectorZ(int, MKL_INT, MKL_Complex16*, MKL_INT, char*);
void PrintArrayS(CBLAS_ORDER*, int, int, MKL_INT*, MKL_INT*, float*, MKL_INT*, char*);
void PrintArrayD(CBLAS_ORDER*, int, int, MKL_INT*, MKL_INT*, double*, MKL_INT*, char*);
void PrintArrayC(CBLAS_ORDER*, int, int, MKL_INT*, MKL_INT*, MKL_Complex8*, MKL_INT*, char*);
void PrintArrayZ(CBLAS_ORDER*, int, int, MKL_INT*, MKL_INT*, MKL_Complex16*, MKL_INT*, char*);
void PrintBandArrayS(CBLAS_ORDER*, int, MKL_INT, MKL_INT, MKL_INT, MKL_INT, float*, MKL_INT, char*);
void PrintBandArrayD(CBLAS_ORDER*, int, MKL_INT, MKL_INT, MKL_INT, MKL_INT, double*, MKL_INT, char*);
void PrintBandArrayC(CBLAS_ORDER*, int, MKL_INT, MKL_INT, MKL_INT, MKL_INT, MKL_Complex8*, MKL_INT, char*);
void PrintBandArrayZ(CBLAS_ORDER*, int, MKL_INT, MKL_INT, MKL_INT, MKL_INT, MKL_Complex16*, MKL_INT, char*);
#ifdef __cplusplus
}
#endif /* __cplusplus */
#endif /* _MKL_EXAMPLE_H_ */
spgemm/mkl_shmem/mklspmm.c → spgemm/mkl_shmem/mklspgemm.c
View file @ 3796eb48
......@@ -5,17 +5,17 @@
#include <mkl.h>
#include "mmio.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';
double time_mic_numeric_mm;
double time_mic_symbolic_mm;
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;
......@@ -26,56 +26,33 @@ typedef struct csr_t {
csi* p;
csi* j;
csv* x;
int numInMats; // number of input matrices stored in this struct
int numRows;
int nnonempty_rows;
int numCols;
int totalnnz;
// int totalnnzmax;
int* nItems; // nnz of each matrix
int* pirow;
int size_pirow;
int* pnr;
int* h_rowDelimiters;
int size_h_rowDelimiters;
int* h_cols;
csv* h_val;
} 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{{{*/
//#include "../../csr/sspmxm.h"
#define CS_MAX(a,b) (((a) > (b)) ? (a) : (b))
csr *csr_spfree(csr *A);
/* free workspace and return a sparse matrix result */
csr *csr_done(csr *C, void *w, void *x, csi ok)
{
//#ifdef KADEBUG
// printf("%d:%s:%s():%d, %s(%lu,%lu,%lu)\n", me, __FILE__, __FUNCTION__,
// __LINE__, C->name, C->m, C->n, C->nzmax);
//#endif
// cs_free(w); /* free workspace */
// cs_free(x);
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)
{
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 *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 */
......@@ -87,8 +64,7 @@ void *csr_realloc(void *p, csi n, size_t size, csi *ok)
}
/* wrapper for realloc */
void *cs_realloc(void *p, csi n, size_t size, csi *ok)
{
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 */
......@@ -98,11 +74,8 @@ void *cs_realloc(void *p, csi n, size_t size, csi *ok)
return((*ok) ? pnew : p); /* return original p if failure */
}
/* change the max # of entries sparse matrix */
csi csr_sprealloc(csr *A, csi nzmax)
{
//printf("%d: nzmax=%d\n", __LINE__, nzmax);
csi csr_sprealloc(csr *A, csi nzmax) {
csi ok, oki = 0, okj = 1, okx = 1;
if (!A)
return(0);
......@@ -111,20 +84,14 @@ csi csr_sprealloc(csr *A, csi nzmax)
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);
//printf("%d: ok=%d oki=%d okj=%d okx=%d\n", __LINE__, ok, oki, okj, okx);
ok = (oki && okj && okx);
//printf("%d: ok=%d oki=%d okj=%d okx=%d\n", __LINE__, ok, oki, okj, okx);
//printf("%s:%d: A->nzmax=%d nzmax=%d\n", __FUNCTION__, __LINE__, A->nzmax, nzmax);
if (ok)
A->nzmax = nzmax;
//printf("%s:%d: A->nzmax=%d nzmax=%d\n", __FUNCTION__, __LINE__, A->nzmax, nzmax);
return(ok);
}
/* free a sparse matrix */
csr *csr_spfree(csr *A)
{
csr *csr_spfree(csr *A) {
if (!A)
return(NULL); /* do nothing if A already NULL */
cs_free(A->p);
......@@ -135,17 +102,13 @@ csr *csr_spfree(csr *A)
A->x = NULL;
cs_free(A->r);
A->r = NULL;
// cs_free(A->name);
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) // floatType f is to suppress compile error
{
//printf("m=%d n=%d nzmax=%d\n", m, n, nzmax);
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 */
// A->name = calloc(MTX_NAME_LEN, sizeof(char));
// sprintf(A->name, "N/A");
if (!A) {
perror("sparse allocation failed");
return(NULL); /* out of memory */
......@@ -153,110 +116,59 @@ csr *csr_spalloc(csi m, csi n, csi nzmax, int values, int triplet, csv f) // flo
A->m = m; /* define dimensions and nzmax */
A->n = n;
A->nzmax = nzmax = CS_MAX(nzmax, 0);
//printf("A m=%d n=%d nzmax=%d\n", A->m, A->n, A->nzmax);
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);
}/*}}}*/
/*csr_multiply{{{*/
/* C = A*B
If nonzero counts of matrices are zero or number of rows/cols is zero, it is advised not to call this function.
*/
// run 8 CUT 0.1 NNZ MTX ~/matrices/nlpkkt200.mtx ~/matrices/nlpkkt200.mtx tmp stdout MNAC SERIAL
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)
{
/** 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;
//printf("%s:%d\n", __FILE__, __LINE__);
if (An != Bm)
return(NULL);
//printf("%s:%d\n", __FILE__, __LINE__);
if (Anzmax == 0 || Bnzmax == 0) {
C = csr_spalloc(Am, Bn, 0, values, 0, tf);
/*#ifdef KADEBUG
printf("%d:%s:%s():%d, A(%s,%lu,%lu,%lu), B(%s,%lu,%lu,%lu), \
One of matrice is Zero! C(%lu,%lu,%lu)\n", me, __FILE__, __FUNCTION__,
__LINE__, A->name, A->m, A->n, A->nzmax, B->name, B->m, B->n,
B->nzmax, C->m, C->n, C->nzmax);
#endif*/
return C;
}
// printf("%s:%d\n", __FILE__, __LINE__);
m = Am;
anz = Ap[Am];
n = Bn;
bnz = Bp[Bm];
for(i = 0; i < n; i++) xb[i] = 0;
//csi* xb = ka_calloc(n, sizeof(csi)); /* get workspace */
for(i = 0; i < n; i++)
xb[i] = 0;
values = (Ax != NULL) && (Bx != NULL);
//csv* x = values ? ka_calloc(n, sizeof(csv)) : NULL; /* get workspace */
csi tnz = (anz + bnz) * 2;
//printf("A->m=%d B->n=%d tnz=%d\n", A->m, B->n, tnz);
C = csr_spalloc(m, n, tnz, values, 0, tf); /* allocate result */
// sprintf(C->name, "C=(%s)*(%s)", A->name, B->name);
/*#ifdef KADEBUG
printf("%d:%s:%s():%d, A(%s,%lu,%lu,%lu), B(%s,%lu,%lu,%lu), \
allocated C(%lu,%lu,%lu)\n", me, __FILE__, __FUNCTION__,
__LINE__, A->name, A->m, A->n, A->nzmax, B->name, B->m, B->n,
B->nzmax, C->m, C->n, C->nzmax);
#endif*/
if (!C || !xb || (values && !x))
return (csr_done(C, xb, x, 0));
// csi zero = 0;
// for(i = 0; i < n; i++)
// xb[i] = zero;
Cp = C->p;
//csi* oldj = C->j;
//printf("C->m=%d C->n=%d C->nzmax=%d C->j[0]=%d\n", C->m, C->n, C->nzmax, C->j[0]);
for (i = 0; i < m; i++) {
//C->j != oldj ? printf("Changed old=%u new=%u\n", oldj, C->j):0;
if ( ( (nz + n) > C->nzmax ) ) {
if(!csr_sprealloc(C, (2 * (C->nzmax) + n) ) ) {
/*#ifdef KADEBUG
printf("%d:%s:%s():%d, A(%s,%lu,%lu,%lu), B(%s,%lu,%lu,%lu), CANNOT BE\
enlargened C(%lu,%lu,%lu)\n", me, __FILE__, __FUNCTION__,
__LINE__, A->name, A->m, A->n, A->nzmax, B->name, B->m, B->n,
B->nzmax, C->m, C->n, C->nzmax);
#endif */
return (csr_done(C, xb, x, 0)); // out of memory
} else {
/*#ifdef KADEBUG
printf("%d:%s:%s():%d, A(%s,%lu,%lu,%lu), B(%s,%lu,%lu,%lu), enlargened \
C(%lu,%lu,%lu)\n", me, __FILE__, __FUNCTION__, __LINE__,
A->name, A->m, A->n, A->nzmax, B->name, B->m, B->n,
B->nzmax, C->m, C->n, C->nzmax);
#endif */
}
//printf("After expansion, C->nzmax=%d\n", C->nzmax);
}
Cj = C->j;
Cx = C->x; /* C->j and C->x may be reallocated */
//printf("i=%d C->j[0]=%d\n", i, C->j[0]);
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(k > n || k < 0)printf("k=%d kp=%d n=%d: %d\n", k, kp, n,__LINE__);
if (xb[k] != i + 1) {
xb[k] = i + 1; /* i is new entry in column j */
// if(nz > C->nzmax) {