vector.c

/*
 * vector.c
 *
 *  routines for handling a distributed vector
 */

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

/* build a row or column distributed vector */
Lalib_vector *
lalib_build_vector(int length, int type, int blocklen, Lalib_grid *grid) {
  Lalib_vector *  vec;
  int             local_length;
  
  vec = (Lalib_vector *) calloc(1, sizeof(Lalib_vector));
  vec->length       = length;
  vec->type         = type;
  vec->grid         = grid;
  vec->blocklength  = blocklen;
  
  /* get number of local vector elements */
  local_length = lalib_vector_local_length(vec);
  
  vec->local_length = local_length;
  vec->data         = (double *) calloc(local_length, sizeof(double));
  
  return(vec);
}

/*
 *  functions describing block scatter distribution of vectors
 */

void lalib_vector_local_index(Lalib_vector *v, int m, int *p, int *i) {
  int   blocklen, nproc, t;
   
  blocklen = v->blocklength;
  nproc = (v->type == LALIB_ROW) ? v->grid->cols : v->grid->rows;
  
  t = nproc*blocklen;
  
  *p = (m % t)/blocklen;
  *i = blocklen * (m / t) + (m % t) % blocklen;
}

void lalib_vector_global_index(Lalib_vector *v, int p, int i, int *m) {
  int   blocklen, nproc, t;
  
  blocklen = v->blocklength;
  nproc = (v->type == LALIB_ROW) ? v->grid->cols : v->grid->rows;
  
  t = nproc*blocklen;
  
  *m = t * (i/blocklen) + blocklen*p + i % blocklen;
}


static int lalib_vector_local_length(Lalib_vector *v) {
  int   blocklen, nproc, len;
  int   t, p, p0;
  
  blocklen = v->blocklength;
  len      = v->length;
  if (v->type == LALIB_ROW) {
    nproc = v->grid->cols;
    p     = v->grid->mycol;
  } else {                 /* v->type == LALIB_COL */
    nproc = v->grid->rows;
    p     = v->grid->myrow;
  }
  
  t = nproc*blocklen;
  p0 = (len % t) / blocklen;
  
  if (p < p0) {
    return(blocklen * ((len/t) + 1));
  } else if (p == p0) {
    return(blocklen * (len/t) + len % blocklen);
  } else {         /* p > p0 */
    return(blocklen * (len/t));
  }
}

/*
 *  compute scalar product of two distributed vectors
 */

static int check_vectors(Lalib_vector *a, Lalib_vector *b);

int lalib_dot_nostride(Lalib_vector *a, Lalib_vector *b, double *result) {
  /* every process sums up its local parts, */
  /* total sum with allreduce across column/row according to vector type */
  
  int       i;
  int       rc;
  double    local_sum = 0;
  MPI_Comm  comm;

  if ( (rc = check_vectors(a, b)) != 0 ) {
    return(rc);
  }
  
  /* sum up local parts */
  for (i = 0; i < a->local_length; i++) {
    local_sum += a->data[i] * b->data[i];
  }

  /* sum up local partial sums */
  comm = (a->type == LALIB_ROW) ? a->grid->row_comm : a->grid->col_comm;
  MPI_Allreduce(&local_sum, result, 1, MPI_DOUBLE, MPI_SUM, comm);
  
  return(0);
}

int lalib_dot_stride(Lalib_vector *a, Lalib_vector *b, double *result) {
  /* every process sums up a strided part of its local data, */
  /* total sum with allreduce across complete grid */

  int       i;
  int       rc;
  double    local_sum = 0;
  MPI_Comm  comm;

  if ( (rc = check_vectors(a, b)) != 0 ) {
    return(rc);
  }

  /* sum up local parts */
  if (a->type == LALIB_ROW) {
    for (i = a->grid->myrow; i < a->local_length; i += a->grid->rows) {
      local_sum += a->data[i] * b->data[i];
    }
  } else {          /* a->type == LALIB_COL */
    for (i = a->grid->mycol; i < a->local_length; i += a->grid->cols) {
      local_sum += a->data[i] * b->data[i];
    }
  }
  
  /* sum up local partial sums */
  MPI_Allreduce(&local_sum, result, 1, MPI_DOUBLE, MPI_SUM, 
		a->grid->grid_comm);
  
  return(0);
}

static int  check_vectors(Lalib_vector *a, Lalib_vector *b) {
  /* check for equal lenghts */
  if (a->length != b->length) {
    printf("unequal vector lengths in lalib_add_vector\n");
    return(-1);
  }

  /* check for compatible types */
  if (a->type != b->type) {
    printf("incompatible vector types in lalib_add_vector\n");
    return(-1);
  }

  /* check for equal blocksize */
  if (a->blocklength != b-> blocklength) {
    printf("unequal block lengths in lalib_add_vector\n");
    return(-1);
  }
   
  /* everything ok */
  return(0);
}