SparseProduct.c

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/// #include "InputOutput.h"
#include "ScalarVectors.h"
#include "hb_io.h"
#include "SparseProduct.h"

/*********************************************************************************/

// This routine creates a sparseMatrix from the next parameters
// * numR defines the number of rows
// * numC defines the number of columns
// * numE defines the number of nonzero elements
// * msr indicates if the MSR is the format used to the sparse matrix
// If msr is actived, numE doesn't include the diagonal elements
// The parameter index indicates if 0-indexing or 1-indexing is used.
void CreateSparseMatrix (ptr_SparseMatrix p_spr, int index, int numR, int numC, int numE, int msr) {
//	printf (" index = %d , numR = %d , numC = %d , numE = %d\n", index, numR, numC, numE);
	// The scalar components of the structure are initiated
	p_spr->dim1 = numR; p_spr->dim2 = numC; 
	// Only one malloc is made for the vectors of indices
	CreateInts (&(p_spr->vptr), numE+numR+1);
	// The first component of the vectors depends on the used format
	*(p_spr->vptr) = ((msr)? (numR+1): 0) + index;
	p_spr->vpos = p_spr->vptr + ((msr)? 0: (numR+1));
	// The number of nonzero elements depends on the format used
	CreateDoubles (&(p_spr->vval), numE+(numR+1)*msr);
}


// This routine creates the first part of a sparseMatrix from the next parameters
// * numR defines the number of rows
// * numC defines the number of columns
// * msr indicates if the MSR is the format used to the sparse matrix
// If msr is actived, numE doesn't include the diagonal elements
// The parameter index indicates if 0-indexing or 1-indexing is used.
void CreateSparseMatrixVptr (ptr_SparseMatrix spr, int numR, int numC, 
												int msr)
	{
		spr->dim1 = numR; spr->dim2 = numC; 
		CreateInts (&(spr->vptr), numR+1);
		*(spr->vptr) = ((msr)? (numR+1): 0);
	}

// This routine creates the second part of a sparseMatrix from the next parameters
// * numR defines the number of rows
// * numC defines the number of columns
// * numE defines the number of nonzero elements
// * msr indicates if the MSR is the format used to the sparse matrix
// If msr is actived, numE doesn't include the diagonal elements
// The parameter index indicates if 0-indexing or 1-indexing is used.
void CreateSparseMatrixValues (ptr_SparseMatrix spr, int numR, int numC, int numE, 
												int msr)
	{
		CreateInts (&(spr->vpos), numE+(numR+1)*msr);
		CreateDoubles (&(spr->vval), numE+(numR+1)*msr);
	}

// This routine liberates the memory related to matrix spr
void RemoveSparseMatrix (ptr_SparseMatrix spr) {
	// First the scalar are initiated
	spr->dim1 = -1; spr->dim2 = -1; 
	// The vectors are liberated
	RemoveInts (&(spr->vptr)); RemoveDoubles (&(spr->vval)); 
}

// This routine liberates the memory related to matrix spr when
// vptr and vpos have been allocated separetely
void RemoveSparseMatrix2 (ptr_SparseMatrix spr)
	{
		spr->dim1 = -1; spr->dim2 = -1; 
		RemoveInts (&(spr->vptr));
		RemoveInts (&(spr->vpos));
		RemoveDoubles (&(spr->vval)); 
	}

/*********************************************************************************/

// This routine creates de sparse matrix dst from the symmetric matrix spr.
// The parameters indexS and indexD indicate, respectivaly, if 0-indexing or 1-indexing is used
// to store the sparse matrices.
void DesymmetrizeSparseMatrices (SparseMatrix src, int indexS, ptr_SparseMatrix dst, int indexD) {
	int n = src.dim1, nnz = 0;
	int *sizes = NULL;
	int *pp1 = NULL, *pp2 = NULL, *pp3 = NULL, *pp4 = NULL, *pp5 = NULL;
	int i, j, dim, indexDS = indexD - indexS;
	double *pd3 = NULL, *pd4 = NULL;

	// The vector sizes is created and initiated
	CreateInts (&sizes, n); InitInts (sizes, n, 0, 0);
	// This loop counts the number of elements in each row
	pp1 = src.vptr; pp3 = src.vpos + *pp1 - indexS;
	pp2 = pp1 + 1 ; pp4 = sizes - indexS;
	for (i=indexS; i<(n+indexS); i++) {
		// The size of the corresponding row is accumulated
		dim = (*pp2 - *pp1); pp4[i] += dim;
		// Now each component of the row is analyzed
		for (j=0; j<dim; j++) {
			// The nondiagonals elements define another element in the graph
			if (*pp3 != i) pp4[*pp3]++;
			pp3++;
		}
		pp1 = pp2++; 
	}
	
	// Compute the number of nonzeros of the new sparse matrix
	nnz = AddInts (sizes, n);
	// Create the new sparse matrix
	CreateSparseMatrix (dst, indexD, n, n, nnz, 0);
	// Fill the vector of pointers
	CopyInts (sizes, (dst->vptr) + 1, n);
	dst->vptr[0] = indexD; TransformLengthtoHeader (dst->vptr, n);
	// The vector sizes is initiated with the beginning of each row
	CopyInts (dst->vptr, sizes, n);
	// This loop fills the contents of vector vpos
	pp1 = src.vptr; pp3 = src.vpos + *pp1 - indexS; 
	pp2 = pp1 + 1 ; pp4 = dst->vpos - indexD; pp5 = sizes - indexS;
	pd3 = src.vval  + *pp1 - indexS; pd4 = dst->vval - indexD;
	for (i=indexS; i<(n+indexS); i++) {
		dim = (*pp2 - *pp1);
		for (j=0; j<dim; j++) {
			// The elements in the i-th row
			pp4[pp5[i]  ] = *pp3+indexDS; 
			pd4[pp5[i]++] = *pd3; 
			if (*pp3 != i) {
				// The nondiagonals elements define another element in the graph
				pp4[pp5[*pp3]  ] = i+indexDS;
				pd4[pp5[*pp3]++] = *pd3;
			}
			pp3++; pd3++;
		}
		pp1 = pp2++;
	}
	// The memory related to the vector sizes is liberated
	RemoveInts (&sizes);
}

/*********************************************************************************/

// This routine creates de sparse matrix dst from the matrix spr.
// The parameters indexS and indexD indicate, respectivaly, if 0-indexing or 1-indexing is used
// to store the sparse matrices.
void TransposeSparseMatrices (SparseMatrix src, int indexS, ptr_SparseMatrix dst, int indexD) {
    int n = src.dim1, nnz = 0;
    int *sizes = NULL;
    int *pp1 = NULL, *pp2 = NULL, *pp3 = NULL, *pp4 = NULL, *pp5 = NULL;
    int i, j, dim, indexDS = indexD - indexS;
    double *pd3 = NULL, *pd4 = NULL;

    // The vector sizes is created and initiated
    CreateInts (&sizes, n); InitInts (sizes, n, 0, 0);
    // This loop counts the number of elements in each row
    pp1 = src.vptr; pp3 = src.vpos + *pp1 - indexS;
    pp2 = pp1 + 1 ; pp4 = sizes - indexS;
    for (i=indexS; i<(n+indexS); i++) {
        // The size of the corresponding row is accumulated
        dim = (*pp2 - *pp1); 
        // Now each component of the row is analyzed
        for (j=0; j<dim; j++) {
            pp4[*pp3]++;
            pp3++;
        }
        pp1 = pp2++; 
    }

    // Compute the number of nonzeros of the new sparse matrix
    nnz = AddInts (sizes, n); 
    // Create the new sparse matrix
    CreateSparseMatrix (dst, indexD, n, n, nnz, 0);
    // Fill the vector of pointers
    CopyInts (sizes, (dst->vptr) + 1, n);
    dst->vptr[0] = indexD; TransformLengthtoHeader (dst->vptr, n);
    // The vector sizes is initiated with the beginning of each row
    CopyInts (dst->vptr, sizes, n);
    // This loop fills the contents of vector vpos
    pp1 = src.vptr; pp3 = src.vpos + *pp1 - indexS; 
    pp2 = pp1 + 1 ; pp4 = dst->vpos - indexD; pp5 = sizes - indexS;
    pd3 = src.vval  + *pp1 - indexS; pd4 = dst->vval - indexD;
    for (i=indexS; i<(n+indexS); i++) {
        dim = (*pp2 - *pp1);
        for (j=0; j<dim; j++) {
            // The elements in the i-th column
            pp4[pp5[*pp3]  ] = i+indexDS;
            pd4[pp5[*pp3]++] = *pd3;
            pp3++; pd3++;
        }
        pp1 = pp2++;
    }
    // The memory related to the vector sizes is liberated
    RemoveInts (&sizes);
}

/*********************************************************************************/

int ReadMatrixHB (char *filename, ptr_SparseMatrix p_spr) {
  int *colptr = NULL;
  double *exact = NULL;
  double *guess = NULL;
  int indcrd;
  char *indfmt = NULL;
  FILE *input;
  char *key = NULL;
  char *mxtype = NULL;
  int ncol;
  int neltvl;
  int nnzero;
  int nrhs;
  int nrhsix;
  int nrow;
  int ptrcrd;
  char *ptrfmt = NULL;
  int rhscrd;
  char *rhsfmt = NULL;
  int *rhsind = NULL;
  int *rhsptr = NULL;
  char *rhstyp = NULL;
  double *rhsval = NULL;
  double *rhsvec = NULL;
  int *rowind = NULL;
  char *title = NULL;
  int totcrd;
  int valcrd;
  char *valfmt = NULL;
  double *values = NULL;

	printf ("\nTEST09\n");
	printf ("  HB_FILE_READ reads all the data in an HB file.\n");
	printf ("  HB_FILE_MODULE is the module that stores the data.\n");

	input = fopen (filename, "r");
	if ( !input ) {
		printf ("\n TEST09 - Warning!\n Error opening the file %s .\n", filename);
		return -1;
	}

	hb_file_read ( input, &title, &key, &totcrd, &ptrcrd, &indcrd,
									&valcrd, &rhscrd, &mxtype, &nrow, &ncol, &nnzero, &neltvl,
									&ptrfmt, &indfmt, &valfmt, &rhsfmt, &rhstyp, &nrhs, &nrhsix,
									&colptr, &rowind, &values, &rhsval, &rhsptr, &rhsind, &rhsvec,
									&guess, &exact );
	fclose (input);

	// Conversion Fortran to C
	CopyShiftInts (colptr, colptr, nrow+1, -1);
	CopyShiftInts (rowind, rowind, nnzero, -1);

	//  Data assignment
	p_spr->dim1 = nrow  ; p_spr->dim2 = ncol  ; 
	p_spr->vptr = colptr; p_spr->vpos = rowind; p_spr->vval = values; 

	//  Memory liberation
	free (exact ); free (guess ); free (indfmt);
	free (key   ); free (mxtype); free (ptrfmt);
	free (rhsfmt); free (rhsind); free (rhsptr);
	free (rhstyp); free (rhsval); free (rhsvec);
	free (title ); free (valfmt);
	
	return 0;
}

/*********************************************************************************/

// This routine computes the product { res += spr * vec }.
// The parameter index indicates if 0-indexing or 1-indexing is used,
void ProdSparseMatrixVector2 (SparseMatrix spr, int index, double *vec, double *res) {
	int i, j;
	int *pp1 = spr.vptr, *pp2 = pp1+1, *pi1 = spr.vpos + *pp1 - index;
	double aux, *pvec = vec - index, *pd2 = res;
	double *pd1 = spr.vval + *pp1 - index;

	// If the MSR format is used, first the diagonal has to be processed
	if (spr.vptr == spr.vpos)
		VvecDoubles (1.0, spr.vval, vec, 1.0, res, spr.dim1);

	for (i=0; i<spr.dim1; i++) {
		// The dot product between the row i and the vector vec is computed
		aux = 0.0;
		for (j=*pp1; j<*pp2; j++)
			aux += *(pd1++) * pvec[*(pi1++)];
//		for (j=spr.vptr[i]; j<spr.vptr[i+1]; j++)
//			aux += spr.vval[j] * pvec[spr.vpos[j]];
		// Accumulate the obtained value on the result
		*(pd2++) += aux; pp1 = pp2++;
	}
}

// This routine computes the product { res += spr * vec }.
// The parameter index indicates if 0-indexing or 1-indexing is used,
void ProdSparseMatrixVectorByRows (SparseMatrix spr, int index, double *vec, double *res) {
	int i, j, dim = spr.dim1;
	int *pp1 = spr.vptr, *pi1 = spr.vpos + *pp1 - index;
	double aux, *pvec = vec + *pp1 - index;
	double *pd1 = spr.vval + *pp1 - index;

	// Process all the rows of the matrix
	for (i=0; i<dim; i++) {
		// The dot product between the row i and the vector vec is computed
		aux = 0.0;
		for (j=pp1[i]; j<pp1[i+1]; j++)
			aux = fma(pd1[j], pvec[pi1[j]], aux);
		// Accumulate the obtained value on the result
		res[i] += aux; 
	}
}

// This routine computes the product { res += spr * vec }.
// The parameter index indicates if 0-indexing or 1-indexing is used,
void ProdSparseMatrixVectorByRows_OMP (SparseMatrix spr, int index, double *vec, double *res) {
	int i, j, dim = spr.dim1;
	int *pp1 = spr.vptr, *pi1 = spr.vpos + *pp1 - index;
	double aux, *pvec = vec + *pp1 - index;
	double *pd1 = spr.vval + *pp1 - index;

	// Process all the rows of the matrix
	#pragma omp parallel for private(j, aux)
	for (i=0; i<dim; i++) {
		// The dot product between the row i and the vector vec is computed
		aux = 0.0;
		for (j=pp1[i]; j<pp1[i+1]; j++)
			aux += pd1[j] * pvec[pi1[j]];
		// Accumulate the obtained value on the result
		res[i] += aux; 
	}
}

/*void ProdSparseMatrixVectorByRows_OMPTasks (SparseMatrix spr, int index, double *vec, double *res, int bm) {
	int i, dim = spr.dim1;

	// Process all the rows of the matrix
	//#pragma omp taskloop grainsize(bm) 
	for ( i=0; i<dim; i+=bm) {
		int cs = dim - i;
		int c = cs < bm ? cs : bm;
//	for (i=0; i<dim; i++) {
	  #pragma omp task depend(inout:res[i:i+c-1]) //shared(c)
		{
//	printf("Task SPMV ---- i: %d, c: %d \n", i, c);
		  int *pp1 = spr.vptr, *pi1 = spr.vpos + *pp1 - index;
	    double aux, *pvec = vec + *pp1 - index;
	    double *pd1 = spr.vval + *pp1 - index;
	  	// The dot product between the row i and the vector vec is computed
		  aux = 0.0;
			for(int idx=i; idx < i+c; idx++){
			//	printf("Task SPMV ---- idx: %d\n", idx);
		  	for (int j=pp1[idx]; j<pp1[idx+1]; j++)
			  	aux += pd1[j] * pvec[pi1[j]];
		  	// Accumulate the obtained value on the result
		  	res[idx] += aux; 
	  	}
		}
	}
}
*/

void ProdSparseMatrixVectorByRows_OMPTasks (SparseMatrix spr, int index, double *vec, double *res, int bm) {
	int i, j, idx, dim = spr.dim1;
	int *pp1 = spr.vptr, *pi1 = spr.vpos + *pp1 - index;
	double aux, *pvec = vec + *pp1 - index;
	double *pd1 = spr.vval + *pp1 - index;

	// Process all the rows of the matrix
	#pragma omp taskloop grainsize(bm) 
	for ( i=0; i<dim; i++ ) {
	  	// The dot product between the row i and the vector vec is computed
		  aux = 0.0;
		  for (j=pp1[i]; j<pp1[i+1]; j++)
			  aux += pd1[j] * pvec[pi1[j]];
		  // Accumulate the obtained value on the result
		  res[i] += aux; 
	}
}


/*********************************************************************************/

// This routine computes the product { res += spr * vec }.
// The parameter index indicates if 0-indexing or 1-indexing is used,
void ProdSparseMatrixVectorByCols (SparseMatrix spr, int index, double *vec, double *res) {
	int i, j, dim = spr.dim1;
	int *pp1 = spr.vptr, *pi1 = spr.vpos + *pp1 - index;
	double aux, *pres = res + *pp1 - index;
	double *pd1 = spr.vval + *pp1 - index;

	// Process all the columns of the matrix
	for (i=0; i<dim; i++) {
		// The result is scaled by the column i and the scalar vec[i]
		aux = vec[i];
		for (j=pp1[i]; j<pp1[i+1]; j++)
			pres[pi1[j]] += pd1[j] * aux;
	}
}

// This routine computes the product { res += spr * vec }.
// The parameter index indicates if 0-indexing or 1-indexing is used,
void ProdSparseMatrixVectorByCols_OMP (SparseMatrix spr, int index, double *vec, double *res) {
	int i, j, dim = spr.dim1;
	int *pp1 = spr.vptr, *pi1 = spr.vpos + *pp1 - index;
	double aux, *pres = res + *pp1 - index;
	double *pd1 = spr.vval + *pp1 - index;

	// Process all the columns of the matrix
	#pragma omp parallel for private(j, aux)
	for (i=0; i<dim; i++) {
		// The result is scaled by the column i and the scalar vec[i]
		for (j=pp1[i]; j<pp1[i+1]; j++) {
			aux = vec[i] * pd1[j];
			#pragma omp atomic
				pres[pi1[j]] += aux;
		}
	}
}

/*********************************************************************************/

void GetDiagonalSparseMatrix2 (SparseMatrix spr, int shft, double *diag, int *posd) {
    int i, j, dim = (spr.dim1 < spr.dim2) ? spr.dim1 : spr.dim2;
    int *pp1 = NULL, *pp2 = NULL, *pi1 = NULL, *pi2 = posd; 
    double *pd1 = NULL, *pd2 = diag;

    if (spr.vptr == spr.vpos)
        CopyDoubles (spr.vval, diag, spr.dim1);
    else {
        pp1 = spr.vptr; pp2 = pp1+1; j = (*pp2-*pp1);
        pi1 = spr.vpos+*pp1; pd1 = spr.vval+*pp1; 
        for (i=0; i<dim; i++) {
            while ((j > 0) && (*pi1 < (i+shft))) {
                pi1++; pd1++; j--;
            }
            *(pd2++) = ((j > 0) && (*pi1 == (i+shft))) ? *pd1: 0.0;
            //*(pi2++) = ((j > 0) && (*pi1 == (i+shft))) ? *pp2-j: -1;
            pi1 += j; pd1 += j; pp1 = (pp2++); j = (*pp2-*pp1);
        }
    }
}

/*********************************************************************************/