/* * 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); }