spmat.c

Go to the documentation of this file.

#include <string.h>
#include "def.h"
#include "blaslapack.h"
#include "struct.h"
#include "internal_proto.h"

void csrcsc(int OUTINDEX, const int nrow, const int ncol, int job,
    double *a, int *ja, int *ia,
    double *ao, int *jao, int *iao) {
  int i,k;
  for (i=0; i<ncol+1; i++) {
    iao[i] = 0;
  }
  // compute nnz of columns of A
  for (i=0; i<nrow; i++) {
    for (k=ia[i]; k<ia[i+1]; k++) {
      iao[ja[k]+1] ++;
    }
  }
  // compute pointers from lengths
  for (i=0; i<ncol; i++) {
    iao[i+1] += iao[i];
  }
  // now do the actual copying
  for (i=0; i<nrow; i++) {
    for (k=ia[i]; k<ia[i+1]; k++) {
      int j = ja[k];
      if (job) {
        ao[iao[j]] = a[k];
      }
      jao[iao[j]++] = i + OUTINDEX;
    }
  }
  /*---- reshift iao and leave */
  for (i=ncol; i>0; i--) {
    iao[i] = iao[i-1] + OUTINDEX;
  }
  iao[0] = OUTINDEX;
}

void sortrow(csrMat *A) {
  /*-------------------------------------------*/
  int nrows = A->nrows;
  int ncols = A->ncols;
  int nnz = A->ia[nrows];
  // work array
  double *b;
  int *jb, *ib;
  Malloc(b, nnz, double);
  Malloc(jb, nnz, int);
  Malloc(ib, ncols+1, int);
  // double transposition
  csrcsc(0, nrows, ncols, 1, A->a, A->ja, A->ia, b, jb, ib);
  csrcsc(0, ncols, nrows, 1, b, jb, ib, A->a, A->ja, A->ia);
  // free
  free(b);
  free(jb);
  free(ib);
}

void csr_resize(int nrow, int ncol, int nnz, csrMat *csr) {
  csr->owndata = 1;
  csr->nrows = nrow;
  csr->ncols = ncol;
  Malloc(csr->ia, nrow+1, int);
  Malloc(csr->ja, nnz, int);
  Malloc(csr->a, nnz, double);
}

void free_csr(csrMat *csr) {
  /* if it does not own the data, do nothing */
  if (!csr->owndata) {
    return;
  }
  free(csr->ia);
  free(csr->ja);
  free(csr->a);
}

void csr_copy(csrMat *A, csrMat *B, int allocB) {
  int nrows = A->nrows;
  int ncols = A->ncols;
  int nnz = A->ia[nrows];
  if (allocB) {
    /* allocate memory */
    csr_resize(nrows, ncols, nnz, B);
  } else {
    /* just set the sizes */
    B->nrows = nrows;
    B->ncols = ncols;
  }
  memcpy(B->ia, A->ia, (nrows+1)*sizeof(int));;
  memcpy(B->ja, A->ja, nnz*sizeof(int));
  memcpy(B->a,  A->a,  nnz*sizeof(double));
}

void free_coo(cooMat *coo) {
  free(coo->ir);
  free(coo->jc);
  free(coo->vv);
}

int cooMat_to_csrMat(int cooidx, cooMat *coo, csrMat *csr) {
  const int nnz = coo->nnz;
  //printf("@@@@ coo2csr, nnz %d\n", nnz);
  /* allocate memory */
  csr_resize(coo->nrows, coo->ncols, nnz, csr);
  const int nrows = coo->nrows;
  /* fill (ia, ja, a) */
  int i;
  for (i=0; i<nrows+1; i++) {
    csr->ia[i] = 0;
  }
  for (i=0; i<nnz; i++) {
    int row = coo->ir[i] - cooidx;
    csr->ia[row+1] ++;
  }
  for (i=0; i<nrows; i++) {
    csr->ia[i+1] += csr->ia[i];
  }
  for (i=0; i<nnz; i++) {
    int row = coo->ir[i] - cooidx;
    int col = coo->jc[i] - cooidx;
    double val = coo->vv[i];
    int k = csr->ia[row];
    csr->a[k] = val;
    csr->ja[k] = col;
    csr->ia[row]++;
  }
  for (i=nrows; i>0; i--) {
    csr->ia[i] = csr->ia[i-1];
  }
  csr->ia[0] = 0;
  /* sort rows ? */
  sortrow(csr);
  return 0;
}


#if 0

double dcsrinfnrm(csrMat *A){
  // computes the inf-norm of A: max abs row sum
  double ta = 0.0;
  int nrows =A->nrows, i, j;
  int *ia = A->ia;
  double *aa = A->a;
  /* for each row */
  for (i=0; i<nrows; i++) {
    /* abs row sum */
    double t = 0.0;
    for (j=ia[i]; j<ia[i+1]; j++) {
      t += fabs(aa[j]);
    }
    /* take max */
    ta = max(ta,t);
  }
  return (ta);
}
#endif

void dcsrmv(char trans, int nrow, int ncol, double *a, 
    int *ia, int *ja, double *x, double *y) {
  int  len, jj=nrow;
  if (trans == 'N') {  
    //#pragma omp parallel for schedule(guided)
    double r;
    /*for (i=0; i<nrow; i++) {
      r = 0.0;
      for (j=ia[i]; j<ia[i+1]; j++) {
      r += a[j] * x[ja[j]];
      }
      y[i] = r; 
      }
      */
    while (jj--) {
      len = *(ia+1) - *ia;
      ia++; 
      r = 0.0;
      while (len--)
        r += x[*ja++]*(*a++);
      *y++ = r;
    }
  } else {
    double xi; 
    int jj, len;
    jj = nrow;
    //-------------------- this is from the matvec used in FILTLAN
    //                     column oriented - gains up to 15% in time 
    double *w = y;
    while (jj--)
      *w++ = 0.0;
    jj = nrow;
    //-------------------- column loop
    while (jj--){
      len = *(ia+1) - *ia;
      ia++;
      xi = *x++;
      while(len--)       
        y[*ja++] += xi*(*a++);
    }
  }
}

void matvec_csr(double *x, double *y, void *data) {
  csrMat *A = (csrMat *) data;
#ifdef USE_MKL
  char cN = 'N';
  /*
  double alp = 1.0, bet = 0.0;
  mkl_dcsrmv(&cN, &(A->nrows), &(A->ncols), &alp, "GXXCXX", 
             A->a, A->ja, A->ia, A->ia+1, x, &bet, y);
  */
  mkl_cspblas_dcsrgemv(&cN, &A->nrows, A->a, A->ia, A->ja, x, y);
#else
  dcsrmv('N', A->nrows, A->ncols, A->a, A->ia, A->ja, x, y);
#endif
}


inline void matadd_insert(double t, csrMat *A, csrMat *C, int i, int *k, 
                          int *j, int *map) {
  /* if this entry already exists in C: 
   * checking if it is the first entry of this row
   * and if column indices match 
   * NOTE that pointer j or k will be modified */
  if (*k > C->ia[i] && C->ja[(*k)-1] == A->ja[*j]) {
    if (map) {
      /* j maps to k-1 in C */
      map[(*j)] = (*k) - 1;
    }
    /* add to existing entry */
    C->a[(*k)-1] += t * A->a[*j];
  } else {
    if (map) {
      /* jA maps to k in C */
      map[*j] = *k; 
    }
    /* create new entry */
    C->ja[*k] = A->ja[*j];
    C->a[*k] = t * A->a[*j];
    (*k)++;
  }
  (*j) ++;
}

int matadd(double alp, double bet, csrMat *A, csrMat *B, csrMat *C,
           int *mapA, int *mapB) {
  int nnzA, nnzB, i, jA, jB, k;
  /* check dimension */
  if (A->nrows != B->nrows || A->ncols != B->ncols) {
    return 1;
  }
  /* nnz of A and B */
  nnzA = A->ia[A->nrows];
  nnzB = B->ia[B->nrows];
  /* alloc C [at most has nnzC = nnzA + nnzB] */
  csr_resize(A->nrows, A->ncols, nnzA+nnzB, C);
  /* nnz counter of C */
  k = 0; 
  C->ia[0] = 0;
  for (i=0; i<A->nrows; i++) {
    /* open row i of A and B */
    /* merging two sorted list */
    for (jA=A->ia[i], jB=B->ia[i]; ; ) {
      if (jA < A->ia[i+1] && jB < B->ia[i+1]) {
        /* will insert the element with smaller col id */
        if (A->ja[jA] <= B->ja[jB]) {
          /* insert jA */
          matadd_insert(alp, A, C, i, &k, &jA, mapA);
        } else {
          /* instert jB */
          matadd_insert(bet, B, C, i, &k, &jB, mapB);
        }
      } else if (jA == A->ia[i+1]) {
        for (; jB < B->ia[i+1]; ) {
          /* instert jB */
          matadd_insert(bet, B, C, i, &k, &jB, mapB);
        }
        break;
      } else {
        for (; jA < A->ia[i+1]; ) {
          /* insert jA */
          matadd_insert(alp, A, C, i, &k, &jA, mapA);
        }
        break;
      }
    }
    C->ia[i+1] = k;
  }
  Realloc(C->ja, k, int);
  Realloc(C->a, k, double);
  return 0;
}

int speye(int n, csrMat *A) {
  int i;
  csr_resize(n, n, n, A);
  for (i=0; i<n; i++) {
    A->ia[i] = A->ja[i] = i;
    A->a[i] = 1.0;
  }
  A->ia[n] = n;
  return 0;
}

/*
 * Diagonal scaling for A such that A = D^{-1}*A*D^{-1}, i.e.,
 * A(i,j) = A(i,j) / (d(i)*d(j)), d = diag(D)
 * @param[in,out] Coo Matrix to scale 
 * @param[in] d: The vector that contains d(i)
 */
void diagScalCoo(cooMat *A, double *d) {
  int i, row, col, nnz = A->nnz;
  /* diagonal scaling for A */
  for (i=0; i<nnz; i++) {
    row = A->ir[i];
    col = A->jc[i];
    A->vv[i] /= d[row] * d[col];
  }
}

/*
 * Diagonal scaling for A such that A = D^{-1}*A*D^{-1}, i.e.,
 * A(i,j) = A(i,j) / (d(i)*d(j)), d = diag(D)
 * @param[in,out] CSR Matrix to scale 
 * @param[in] d: The vector that contains d(i)
 */
void diagScalCsr(csrMat *A, double *d) {
  int i, j;
  /* diagonal scaling for A */
  for (i=0; i<A->nrows; i++) {
    for (j=A->ia[i]; j<A->ia[i+1]; j++) {
      A->a[j] /= d[i] * d[A->ja[j]];
    }
  }
}

/*
 * Extract the diagonal entries of csr matrix B
 *
 * @param[in]  B Matrix to extract the diagonal
 * @param[out] D preallocated vector of lengeth B.nrows
 */
void extrDiagCsr(csrMat *B, double *d) {
  int i, j, nrows = B->nrows;
  for (i=0; i<nrows; i++) {
    d[i] = 0.0;
    for (j=B->ia[i]; j<B->ia[i+1]; j++) {
      if (i == B->ja[j]) {
        d[i] = B->a[j];
      }
    }
  }
}

/*
 * Extract the upper triangular part of csr matrix A
 *
 * @param[in]  A Matrix 
 * @param[out] U Matrix
 */
void triuCsr(csrMat *A, csrMat *U) {
  EVSL_Int i, j;
  EVSL_Unsigned k = 0, nnzu = 0;

  /* count nnz in U */
  for (i = 0; i < A->nrows; i++) {
    for (j = A->ia[i]; j < A->ia[i+1]; j++) {
      if (i <= A->ja[j]) {
        nnzu++;
      }
    }
  }
  /* allocate memory for U */
  csr_resize(A->nrows, A->ncols, nnzu, U);
  U->ia[0] = 0;
  /* copy entries from A to U */
  for (i = 0; i < A->nrows; i++) {
    for (j = A->ia[i]; j < A->ia[i+1]; j++) {
      if (i <= A->ja[j]) {
        U->ja[k] = A->ja[j];
        U->a[k++] = A->a[j];
      }
    }
    U->ia[i+1] = k;
  }
  CHKERR(k != nnzu);
}