[cig-commits] r6666 - cs/cigma/trunk/sandbox/c

[luis at geodynamics.org]

Author: luis
Date: 2007-04-24 14:59:03 -0700 (Tue, 24 Apr 2007)
New Revision: 6666

Added:
   cs/cigma/trunk/sandbox/c/qr.c
   cs/cigma/trunk/sandbox/c/qr.h
Log:
Provide three ways to specify the quadrature rules (from a text file,
from an hdf5 file, or directly from an array)


Added: cs/cigma/trunk/sandbox/c/qr.c
===================================================================
--- cs/cigma/trunk/sandbox/c/qr.c	2007-04-24 21:52:29 UTC (rev 6665)
+++ cs/cigma/trunk/sandbox/c/qr.c	2007-04-24 21:59:03 UTC (rev 6666)
@@ -0,0 +1,189 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include "qr.h"
+#include "hdf5.h"
+
+
+static int qr_init1_read_weights(rule_t *qr, hid_t dset_id);
+static int qr_init1_read_points(rule_t *qr, hid_t dset_id);
+
+
+int qr_fini(rule_t *qr)
+{
+    free(qr->weights);
+    free(qr->points);
+    return 0;
+}
+
+int qr_init0(rule_t *qr, char *filename)
+{
+    int q, nq;
+    int i, n;
+
+    FILE *fp;
+    int ret;
+    
+    fp = fopen(filename, "r");
+    if (fp == NULL)
+    {
+        return -1;
+    }
+
+    ret = fscanf(fp, "%d", &nq);
+    ret = fscanf(fp, "%d", &n);
+
+    qr->nq = nq;
+    qr->n  = n;
+    qr->weights = (double *)malloc(nq * sizeof(double));
+    qr->points  = (double *)malloc(nq * n * sizeof(double));
+
+    for (q = 0; q < nq; q++)
+    {
+        ret = fscanf(fp, "%lf", &(qr->weights[q]));
+        for (i = 0; i < n; i++)
+        {
+            ret = fscanf(fp, "%lf", &(qr->points[n*q + i]));
+        }
+    }
+
+    fclose(fp);
+    return 0;
+}
+
+int qr_init1(rule_t *qr, char *filename, char *path)
+{
+    hid_t file_id;
+    hid_t group_id;
+    hid_t weights_dset;
+    hid_t points_dset;
+
+    herr_t status;
+    int ierr;
+
+    /* for now, zero out n & nq */
+    qr->nq = 0;
+    qr->n = 0;
+
+    /* open the file in read-only mode */
+    file_id = H5Fopen(filename, H5F_ACC_RDONLY, H5P_DEFAULT);
+    if (file_id < 0)
+    {
+        return -1;
+    }
+
+    /* open group at specified path */
+    group_id = H5Gopen(file_id, path);
+    if (group_id < 0)
+    {
+        H5Fclose(file_id);
+        return -2;
+    }
+
+    /* read quadrature weights */
+    weights_dset = H5Dopen(group_id, "weights");
+    ierr = qr_init1_read_weights(qr, weights_dset);
+    status = H5Dclose(weights_dset);
+
+    /* read quadrature points */
+    points_dset = H5Dopen(group_id, "points");
+    ierr = qr_init1_read_points(qr, points_dset);
+    status = H5Dclose(points_dset);
+
+    /* finalize */
+    status = H5Gclose(group_id);
+    status = H5Fclose(file_id);
+
+    return 0;
+}
+
+
+void qr_apply(rule_t *qr, double (*f)(double *), double *sum)
+{
+    int q;              // for looping over quadrature points
+    int nq = qr->nq;    // number of quadrature points
+    int n  = qr->n;     // dimension of points
+    
+    double *x = qr->points;
+    double *w = qr->weights;
+    double acc = 0.0;
+
+    for (q = 0; q < nq; q++)
+        acc += w[q] * f(&x[q*n]);
+
+    *sum = acc;
+}
+
+
+
+/*
+ * Private functions
+ */
+
+static int qr_init1_read_weights(rule_t *qr, hid_t dset_id)
+{
+    herr_t status;
+    hid_t weights_dspace;
+    int rank;
+
+    weights_dspace = H5Dget_space(dset_id);
+
+    rank = H5Sget_simple_extent_ndims(weights_dspace);
+    if (rank != 1)
+    {
+        return -1;  // error: expecting 1d array
+    }
+
+    if (qr->nq == 0)
+    {
+        qr->nq = H5Sget_simple_extent_npoints(weights_dspace);
+    }
+
+    status = H5Sclose(weights_dspace);
+
+    qr->weights = (double *)malloc((qr->nq) * sizeof(double));
+
+    status = H5Dread(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
+                     H5P_DEFAULT, qr->weights);
+
+    return 0;
+}
+
+static int qr_init1_read_points(rule_t *qr, hid_t dset_id)
+{
+    herr_t status;
+    hid_t points_dspace;
+    hsize_t npoints;
+    hsize_t *dims;
+    int rank;
+
+    points_dspace = H5Dget_space(dset_id);
+
+    rank = H5Sget_simple_extent_ndims(points_dspace);
+    if (rank != 2 && rank != 1)
+    {
+        return -1;
+    }
+
+    dims = (hsize_t *)malloc(rank * sizeof(hsize_t));
+    status = H5Sget_simple_extent_dims(points_dspace, dims, NULL);
+    npoints = H5Sget_simple_extent_npoints(points_dspace);
+    if (qr->nq != 0)
+    {
+        if (qr->nq != dims[0])
+            return -2;
+    }
+    else
+    {
+        qr->nq = dims[0];
+        qr->n  = (rank > 1) ? dims[1] : 1;
+    }
+
+    qr->points = (double *)malloc(npoints * sizeof(double));
+
+    status = H5Dread(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
+                     H5P_DEFAULT, qr->points);
+
+    free(dims);
+
+    return 0;
+}

Added: cs/cigma/trunk/sandbox/c/qr.h
===================================================================
--- cs/cigma/trunk/sandbox/c/qr.h	2007-04-24 21:52:29 UTC (rev 6665)
+++ cs/cigma/trunk/sandbox/c/qr.h	2007-04-24 21:59:03 UTC (rev 6666)
@@ -0,0 +1,24 @@
+#ifndef __CIGMA_QR_H__
+#define __CIGMA_QR_H__
+
+typedef struct
+{
+    int n;
+    int nq;
+    double *weights;
+    double *points;
+} rule_t;
+
+typedef void (*generator_t)(int, int, double *, double *);
+typedef double (*functor_t)(double *);
+
+int qr_init0(rule_t *qr, char *filename);
+int qr_init1(rule_t *qr, char *filename, char *path);
+int qr_init2(rule_t *qr, generator_t gen);
+int qr_fini(rule_t *qr);
+
+void qr_get(rule_t *qr, double *w, double *x);
+
+void qr_apply(rule_t *qr, functor_t f, double *sum);
+
+#endif