| 1 | #include <mpi.h>
|
|---|
| 2 | #ifdef _CIVL
|
|---|
| 3 | #include <stdlib.h>
|
|---|
| 4 | #include <civlc.cvh>
|
|---|
| 5 | #endif
|
|---|
| 6 |
|
|---|
| 7 | #define HYPRE_BigInt int
|
|---|
| 8 |
|
|---|
| 9 | // seq_mv.h :
|
|---|
| 10 |
|
|---|
| 11 | typedef struct
|
|---|
| 12 | {
|
|---|
| 13 | double *data;
|
|---|
| 14 | int size;
|
|---|
| 15 |
|
|---|
| 16 | /* Does the Vector create/destroy `data'? */
|
|---|
| 17 | int owns_data;
|
|---|
| 18 |
|
|---|
| 19 | /* For multivectors...*/
|
|---|
| 20 | int num_vectors; /* the above "size" is size of one vector */
|
|---|
| 21 | int multivec_storage_method;
|
|---|
| 22 | /* ...if 0, store colwise v0[0], v0[1], ..., v1[0], v1[1], ... v2[0]... */
|
|---|
| 23 | /* ...if 1, store rowwise v0[0], v1[0], ..., v0[1], v1[1], ... */
|
|---|
| 24 | /* With colwise storage, vj[i] = data[ j*size + i]
|
|---|
| 25 | With rowwise storage, vj[i] = data[ j + num_vectors*i] */
|
|---|
| 26 | int vecstride, idxstride;
|
|---|
| 27 | /* ... so vj[i] = data[ j*vecstride + i*idxstride ] regardless of row_storage.*/
|
|---|
| 28 |
|
|---|
| 29 | } hypre_Vector;
|
|---|
| 30 |
|
|---|
| 31 | #define hypre_VectorData(vector) ((vector) -> data)
|
|---|
| 32 | #define hypre_VectorSize(vector) ((vector) -> size)
|
|---|
| 33 | #define hypre_VectorOwnsData(vector) ((vector) -> owns_data)
|
|---|
| 34 | #define hypre_VectorNumVectors(vector) ((vector) -> num_vectors)
|
|---|
| 35 | #define hypre_VectorMultiVecStorageMethod(vector) ((vector) -> multivec_storage_method)
|
|---|
| 36 | #define hypre_VectorVectorStride(vector) ((vector) -> vecstride )
|
|---|
| 37 | #define hypre_VectorIndexStride(vector) ((vector) -> idxstride )
|
|---|
| 38 |
|
|---|
| 39 |
|
|---|
| 40 | // vector.c :
|
|---|
| 41 |
|
|---|
| 42 | double hypre_SeqVectorInnerProd( hypre_Vector *x,
|
|---|
| 43 | hypre_Vector *y )
|
|---|
| 44 | {
|
|---|
| 45 | double *x_data = hypre_VectorData(x);
|
|---|
| 46 | double *y_data = hypre_VectorData(y);
|
|---|
| 47 | int size = hypre_VectorSize(x);
|
|---|
| 48 |
|
|---|
| 49 | int i;
|
|---|
| 50 |
|
|---|
| 51 | double result = 0.0;
|
|---|
| 52 |
|
|---|
| 53 | size *=hypre_VectorNumVectors(x);
|
|---|
| 54 |
|
|---|
| 55 | #ifdef HYPRE_USING_OPENMP
|
|---|
| 56 | #pragma omp parallel for private(i) reduction(+:result) schedule(static)
|
|---|
| 57 | #endif
|
|---|
| 58 | for (i = 0; i < size; i++)
|
|---|
| 59 | result += y_data[i] * x_data[i];
|
|---|
| 60 |
|
|---|
| 61 | return result;
|
|---|
| 62 | }
|
|---|
| 63 |
|
|---|
| 64 |
|
|---|
| 65 | // parcsr_mv.h:
|
|---|
| 66 |
|
|---|
| 67 | typedef struct
|
|---|
| 68 | {
|
|---|
| 69 | int length;
|
|---|
| 70 | HYPRE_BigInt row_start;
|
|---|
| 71 | HYPRE_BigInt row_end;
|
|---|
| 72 | int storage_length;
|
|---|
| 73 | int *proc_list;
|
|---|
| 74 | HYPRE_BigInt *row_start_list;
|
|---|
| 75 | HYPRE_BigInt *row_end_list;
|
|---|
| 76 | int *sort_index;
|
|---|
| 77 | } hypre_IJAssumedPart;
|
|---|
| 78 |
|
|---|
| 79 |
|
|---|
| 80 | typedef struct
|
|---|
| 81 | {
|
|---|
| 82 | MPI_Comm comm;
|
|---|
| 83 |
|
|---|
| 84 | HYPRE_BigInt global_size;
|
|---|
| 85 | HYPRE_BigInt first_index;
|
|---|
| 86 | HYPRE_BigInt last_index;
|
|---|
| 87 | HYPRE_BigInt *partitioning;
|
|---|
| 88 | hypre_Vector *local_vector;
|
|---|
| 89 |
|
|---|
| 90 | int owns_data;
|
|---|
| 91 | int owns_partitioning;
|
|---|
| 92 |
|
|---|
| 93 | hypre_IJAssumedPart *assumed_partition;
|
|---|
| 94 | } hypre_ParVector;
|
|---|
| 95 |
|
|---|
| 96 |
|
|---|
| 97 | #define hypre_ParVectorComm(vector) ((vector) -> comm)
|
|---|
| 98 | #define hypre_ParVectorGlobalSize(vector) ((vector) -> global_size)
|
|---|
| 99 | #define hypre_ParVectorFirstIndex(vector) ((vector) -> first_index)
|
|---|
| 100 | #define hypre_ParVectorLastIndex(vector) ((vector) -> last_index)
|
|---|
| 101 | #define hypre_ParVectorPartitioning(vector) ((vector) -> partitioning)
|
|---|
| 102 | #define hypre_ParVectorLocalVector(vector) ((vector) -> local_vector)
|
|---|
| 103 | #define hypre_ParVectorOwnsData(vector) ((vector) -> owns_data)
|
|---|
| 104 | #define hypre_ParVectorOwnsPartitioning(vector) ((vector) -> owns_partitioning)
|
|---|
| 105 | #define hypre_ParVectorNumVectors(vector)\
|
|---|
| 106 | (hypre_VectorNumVectors( hypre_ParVectorLocalVector(vector) ))
|
|---|
| 107 |
|
|---|
| 108 | #define hypre_ParVectorAssumedPartition(vector) ((vector) -> assumed_partition)
|
|---|
| 109 |
|
|---|
| 110 |
|
|---|
| 111 | // par_vector.c :
|
|---|
| 112 |
|
|---|
| 113 | double
|
|---|
| 114 | hypre_ParVectorInnerProd( hypre_ParVector *x,
|
|---|
| 115 | hypre_ParVector *y )
|
|---|
| 116 | {
|
|---|
| 117 | MPI_Comm comm = hypre_ParVectorComm(x);
|
|---|
| 118 | hypre_Vector *x_local = hypre_ParVectorLocalVector(x);
|
|---|
| 119 | hypre_Vector *y_local = hypre_ParVectorLocalVector(y);
|
|---|
| 120 |
|
|---|
| 121 | double result = 0.0;
|
|---|
| 122 | double local_result = hypre_SeqVectorInnerProd(x_local, y_local);
|
|---|
| 123 |
|
|---|
| 124 | MPI_Allreduce(&local_result, &result, 1, MPI_DOUBLE, MPI_SUM, comm);
|
|---|
| 125 |
|
|---|
| 126 | return result;
|
|---|
| 127 | }
|
|---|
| 128 |
|
|---|
| 129 |
|
|---|
| 130 | /* Stripped down driver for AMG2013 parallel inner product routine. */
|
|---|
| 131 | #define XVET x.local_vector
|
|---|
| 132 | #define YVET y.local_vector
|
|---|
| 133 |
|
|---|
| 134 | #ifdef _CIVL
|
|---|
| 135 | $input int VECTOR_LENGTH;
|
|---|
| 136 | $assume(0 <= VECTOR_LENGTH && VECTOR_LENGTH < 10);
|
|---|
| 137 | #endif
|
|---|
| 138 | int main() {
|
|---|
| 139 | hypre_ParVector x, y;
|
|---|
| 140 | int nprocs;
|
|---|
| 141 |
|
|---|
| 142 | MPI_Init(NULL, NULL);
|
|---|
| 143 | MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
|
|---|
| 144 | #ifdef _CIVL
|
|---|
| 145 | x.comm = MPI_COMM_WORLD;
|
|---|
| 146 | y.comm = MPI_COMM_WORLD;
|
|---|
| 147 | XVET = (hypre_Vector *)malloc(sizeof(hypre_Vector));
|
|---|
| 148 | YVET = (hypre_Vector *)malloc(sizeof(hypre_Vector));
|
|---|
| 149 | XVET->data = (double *)malloc(sizeof(double) * VECTOR_LENGTH * nprocs);
|
|---|
| 150 | YVET->data = (double *)malloc(sizeof(double) * VECTOR_LENGTH * nprocs);
|
|---|
| 151 | XVET->size = VECTOR_LENGTH;
|
|---|
| 152 | YVET->size = VECTOR_LENGTH;
|
|---|
| 153 | XVET->num_vectors = nprocs;
|
|---|
| 154 | YVET->num_vectors = nprocs;
|
|---|
| 155 | #endif
|
|---|
| 156 | double result = hypre_ParVectorInnerProd(&x, &y);
|
|---|
| 157 |
|
|---|
| 158 | MPI_Finalize();
|
|---|
| 159 | free(XVET->data);
|
|---|
| 160 | free(YVET->data);
|
|---|
| 161 | free(XVET);
|
|---|
| 162 | free(YVET);
|
|---|
| 163 | #ifdef DEBUG
|
|---|
| 164 | #include <stdio.h>
|
|---|
| 165 | printf("result = %f\n", result);
|
|---|
| 166 | #endif
|
|---|
| 167 | return result != 0;
|
|---|
| 168 | }
|
|---|