Ignore:
Timestamp:
Nov 17, 2004, 11:40:28 PM (20 years ago)
Author:
steve
Message:

testing sparse

File:
1 edited

Legend:

Unmodified
Added
Removed

  • inundation/ga/storm_surge/pyvolution/sparse_ext.c

    r586 r587  
    1414#include "math.h"
    1515
    16 //Shared code snippets
    17 #include "util_ext.h"
    18 #include "quantity_ext.h"
    19 
    20 
    21 
    22 int _compute_csr_mult(int M,
    23                       double* data,
    24                       int* colind,
    25                       int* row_ptr,
    26                       double* R) {
    27 
    28                        
     16int _csr_mv(int M,
     17            double* data,
     18            int* colind,
     19            int* row_ptr,
     20            double* x,
     21            double* y) {
     22               
    2923  int i, j, ckey;
    3024
    31   for (i=0; i<M; ):
    32                 for ckey in range(self.row_ptr[i],self.row_ptr[i+1]):
    33                     j = self.colind[ckey]
    34                     R[i] += self.data[ckey]*B[j]             
     25  for (i=0; i<M; i++ )
     26    for (ckey=row_ptr[i]; ckey<row_ptr[i+1]; ckey++) {
     27      j = colind[ckey];
     28      y[i] += data[ckey]*x[j];
     29    }             
    3530 
    36   for (k=0; k<N; k++) {
    37     index3 = 3*k;
    38    
    39     if (number_of_boundaries[k] < 2) {
    40       //Two or three true neighbours
    41 
    42       //Get indices of neighbours (or self when used as surrogate)         
    43       //k0, k1, k2 = surrogate_neighbours[k,:]
    44      
    45       k0 = surrogate_neighbours[index3 + 0];
    46       k1 = surrogate_neighbours[index3 + 1];
    47       k2 = surrogate_neighbours[index3 + 2];           
    48       if (k0 == k1 || k1 == k2) return -1;     
    49 
    50       //Get data
    51       q0 = centroid_values[k0];
    52       q1 = centroid_values[k1];
    53       q2 = centroid_values[k2];                   
    54 
    55       x0 = centroids[k0*2]; y0 = centroids[k0*2+1];
    56       x1 = centroids[k1*2]; y1 = centroids[k1*2+1];
    57       x2 = centroids[k2*2]; y2 = centroids[k2*2+1];             
    58 
    59       //Gradient
    60       _gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2, &a[k], &b[k]);
    61      
    62     } else if (number_of_boundaries[k] == 2) {
    63       //One true neighbour
    64 
    65       //#Get index of the one neighbour
    66       i=0; k0 = k;
    67       while (i<3 && k0==k) {
    68         k0 = surrogate_neighbours[index3 + i];
    69         i++;
    70       }
    71       if (k0 == k) return -1;
    72      
    73       k1 = k; //self
    74 
    75       //Get data
    76       q0 = centroid_values[k0];
    77       q1 = centroid_values[k1];
    78            
    79       x0 = centroids[k0*2]; y0 = centroids[k0*2+1];
    80       x1 = centroids[k1*2]; y1 = centroids[k1*2+1];
    81 
    82       //Gradient
    83       det = x0*y1 - x1*y0;
    84       if (det != 0.0) {
    85         a[k] = (y1*q0 - y0*q1)/det;
    86         b[k] = (x0*q1 - x1*q0)/det;
    87       }
    88     }
    89     //    else:
    90     //        #No true neighbours -     
    91     //        #Fall back to first order scheme
    92     //        pass
    93   }
    9431  return 0;
    9532}           
    9633
    9734
    98 int _extrapolate(int N,
    99                  double* centroids,
    100                  double* centroid_values,               
    101                  double* vertex_coordinates,
    102                  double* vertex_values,
    103                  double* a,
    104                  double* b) {
     35                     
    10536 
    106   int k, k2, k3, k6;
    107   double x, y, x0, y0, x1, y1, x2, y2;
    108                      
    109   for (k=0; k<N; k++) {
    110     k6 = 6*k;
    111     k3 = 3*k;   
    112     k2 = 2*k;
    113        
    114     //Centroid coordinates           
    115     x = centroids[k2]; y = centroids[k2+1];
     37/////////////////////////////////////////////////
     38// Gateways to Python
     39PyObject *csr_mv(PyObject *self, PyObject *args) {
     40 
     41  PyObject *csr_sparse, // input sparse matrix
     42    *R;                 // output wrapped vector
     43 
     44  PyArrayObject
     45    *data,            //Non Zeros Data array
     46    *colind,          //Column indices array
     47    *row_ptr,         //Row pointers array
     48    *x,               //Input vector array
     49    *y;               //Return vector array
     50   
     51  int dimensions[1], M, err;
     52 
     53  // Convert Python arguments to C 
     54  if (!PyArg_ParseTuple(args, "OO", &csr_sparse, &x))
     55    return NULL;
    11656
    117     //vertex coordinates
    118     //x0, y0, x1, y1, x2, y2 = X[k,:]
    119     x0 = vertex_coordinates[k6 + 0];
    120     y0 = vertex_coordinates[k6 + 1];   
    121     x1 = vertex_coordinates[k6 + 2];
    122     y1 = vertex_coordinates[k6 + 3];       
    123     x2 = vertex_coordinates[k6 + 4];
    124     y2 = vertex_coordinates[k6 + 5];           
    12557
    126     //Extrapolate
    127     vertex_values[k3+0] = centroid_values[k] + a[k]*(x0-x) + b[k]*(y0-y);
    128     vertex_values[k3+1] = centroid_values[k] + a[k]*(x1-x) + b[k]*(y1-y);
    129     vertex_values[k3+2] = centroid_values[k] + a[k]*(x2-x) + b[k]*(y2-y);
     58  data =  (PyArrayObject*)
     59    PyObject_GetAttrString(csr_sparse, "data");     
     60  if (!data)
     61    return NULL; 
     62
     63  colind = (PyArrayObject*)
     64    PyObject_GetAttrString(csr_sparse, "colind");
     65  if (!colind) return NULL;   
     66 
     67  row_ptr = (PyArrayObject*)
     68    PyObject_GetAttrString(csr_sparse, "row_ptr");   
     69  if (!row_ptr) return NULL;       
     70 
     71  M = (row_ptr -> dimensions[0])-1;
    13072   
     73  //Allocate space for return vectors y (don't DECREF)
     74  dimensions[0] = M;
     75  y = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
     76 
     77  err = _csr_mv(M,
     78                (double*) data -> data,
     79                (int*)    colind -> data,
     80                (int*)    row_ptr -> data,
     81                (double*) x -> data,
     82                (double*) y -> data);     
     83                           
     84  if (err != 0) {
     85    PyErr_SetString(PyExc_RuntimeError, "matrix vector mult could not be calculated");
     86    return NULL;
    13187  }
    132   return 0;
     88                     
     89  //Release                 
     90  Py_DECREF(data);   
     91  Py_DECREF(colind);   
     92  Py_DECREF(row_ptr);
     93  Py_DECREF(x);           
     94         
     95  //Build result, release and return
     96  R = Py_BuildValue("O", PyArray_Return(y));
     97  Py_DECREF(y);       
     98  return R;
    13399}
    134100
    … …  
    136102
    137103
    138 int _interpolate(int N,
    139                  double* vertex_values,
    140                  double* edge_values) {
    141  
    142   int k, k3;
    143   double q0, q1, q2;
    144 
    145                      
    146   for (k=0; k<N; k++) {
    147     k3 = 3*k;   
    148    
    149     q0 = vertex_values[k3 + 0];
    150     q1 = vertex_values[k3 + 1];
    151     q2 = vertex_values[k3 + 2];
    152            
    153     //printf("%f, %f, %f\n", q0, q1, q2);
    154     edge_values[k3 + 0] = 0.5*(q1+q2);
    155     edge_values[k3 + 1] = 0.5*(q0+q2); 
    156     edge_values[k3 + 2] = 0.5*(q0+q1);
    157   } 
    158   return 0;
    159 }
    160 
    161 int _update(int N,
    162             double timestep,
    163             double* centroid_values,
    164             double* explicit_update,
    165             double* semi_implicit_update) {
    166   //Update centroid values based on values stored in
    167   //explicit_update and semi_implicit_update as well as given timestep
    168  
    169  
    170   int k;
    171   double denominator, x;
    172 
    173  
    174   //Divide semi_implicit update by conserved quantity 
    175   for (k=0; k<N; k++) { 
    176     x = centroid_values[k];
    177     if (x == 0.0) {
    178       //FIXME: Is this right
    179       semi_implicit_update[k] = 0.0;
    180     } else {
    181       semi_implicit_update[k] /= x;
    182     }   
    183   }             
    184  
    185  
    186   //Explicit updates
    187   for (k=0; k<N; k++) {
    188     centroid_values[k] += timestep*explicit_update[k];
    189   }
    190            
    191   //Semi implicit updates
    192   for (k=0; k<N; k++) { 
    193     denominator = 1.0 - timestep*semi_implicit_update[k]; 
    194    
    195     if (denominator == 0.0) {
    196       return -1;
    197     } else {
    198       //Update conserved_quantities from semi implicit updates
    199       centroid_values[k] /= denominator;
    200     }
    201   }
    202   return 0;
    203 }
    204                      
    205  
    206 /////////////////////////////////////////////////
    207 // Gateways to Python
    208 PyObject *update(PyObject *self, PyObject *args) {
    209  
    210   PyObject *quantity;
    211   PyArrayObject *centroid_values, *explicit_update, *semi_implicit_update;
    212  
    213   double timestep; 
    214   int N, err;
    215  
    216   // Convert Python arguments to C 
    217   if (!PyArg_ParseTuple(args, "Od", &quantity, &timestep))
    218     return NULL;
    219 
    220   centroid_values = (PyArrayObject*)
    221     PyObject_GetAttrString(quantity, "centroid_values");           
    222   if (!centroid_values) return NULL;               
    223 
    224   explicit_update = (PyArrayObject*)
    225     PyObject_GetAttrString(quantity, "explicit_update");   
    226   if (!explicit_update) return NULL;         
    227  
    228   semi_implicit_update = (PyArrayObject*)
    229     PyObject_GetAttrString(quantity, "semi_implicit_update");   
    230   if (!semi_implicit_update) return NULL;           
    231 
    232   N = centroid_values -> dimensions[0]; 
    233  
    234  
    235   err = _update(N, timestep,
    236                 (double*) centroid_values -> data,
    237                 (double*) explicit_update -> data,
    238                 (double*) semi_implicit_update -> data);               
    239                
    240                      
    241   if (err != 0) {
    242     PyErr_SetString(PyExc_RuntimeError,
    243                     "Zero division in semi implicit update - call Stephen :)");
    244     return NULL;
    245   }                 
    246  
    247   //Release and return               
    248   Py_DECREF(centroid_values);   
    249   Py_DECREF(explicit_update);
    250   Py_DECREF(semi_implicit_update);   
    251 
    252   return Py_BuildValue("");
    253 }
    254 
    255 
    256 PyObject *interpolate_from_vertices_to_edges(PyObject *self, PyObject *args) {
    257  
    258   PyObject *quantity;
    259   PyArrayObject *vertex_values, *edge_values;
    260    
    261   int N, err;
    262  
    263   // Convert Python arguments to C 
    264   if (!PyArg_ParseTuple(args, "O", &quantity))
    265     return NULL;
    266 
    267   vertex_values = (PyArrayObject*)
    268     PyObject_GetAttrString(quantity, "vertex_values");           
    269   if (!vertex_values) return NULL;               
    270 
    271   edge_values = (PyArrayObject*)
    272     PyObject_GetAttrString(quantity, "edge_values");   
    273   if (!edge_values) return NULL;         
    274 
    275   N = vertex_values -> dimensions[0]; 
    276  
    277   err = _interpolate(N,
    278                      (double*) vertex_values -> data,
    279                      (double*) edge_values -> data);
    280                      
    281   if (err != 0) {
    282     PyErr_SetString(PyExc_RuntimeError, "Interpolate could not be computed");
    283     return NULL;
    284   }                 
    285  
    286   //Release and return               
    287   Py_DECREF(vertex_values);   
    288   Py_DECREF(edge_values);
    289 
    290   return Py_BuildValue("");
    291 }
    292 
    293 
    294 PyObject *compute_gradients(PyObject *self, PyObject *args) {
    295  
    296   PyObject *quantity, *domain, *R;
    297   PyArrayObject
    298     *centroids,            //Coordinates at centroids
    299     *centroid_values,      //Values at centroids   
    300     *number_of_boundaries, //Number of boundaries for each triangle
    301     *surrogate_neighbours, //True neighbours or - if one missing - self
    302     *a, *b;                //Return values
    303    
    304   int dimensions[1], N, err;
    305  
    306   // Convert Python arguments to C 
    307   if (!PyArg_ParseTuple(args, "O", &quantity))
    308     return NULL;
    309 
    310   domain = PyObject_GetAttrString(quantity, "domain");     
    311   if (!domain)
    312     return NULL; 
    313 
    314   //Get pertinent variables
    315   centroids = (PyArrayObject*)
    316     PyObject_GetAttrString(domain, "centroid_coordinates");
    317   if (!centroids) return NULL;   
    318  
    319   centroid_values = (PyArrayObject*)
    320     PyObject_GetAttrString(quantity, "centroid_values");   
    321   if (!centroid_values) return NULL;       
    322  
    323   surrogate_neighbours = (PyArrayObject*)
    324     PyObject_GetAttrString(domain, "surrogate_neighbours");
    325   if (!surrogate_neighbours) return NULL;       
    326  
    327   number_of_boundaries = (PyArrayObject*)
    328     PyObject_GetAttrString(domain, "number_of_boundaries");           
    329   if (!number_of_boundaries) return NULL;           
    330  
    331   N = centroid_values -> dimensions[0];
    332    
    333   //Release
    334   Py_DECREF(domain);     
    335          
    336   //Allocate space for return vectors a and b (don't DECREF)
    337   dimensions[0] = N;
    338   a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
    339   b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); 
    340        
    341 
    342  
    343   err = _compute_gradients(N,
    344                            (double*) centroids -> data,
    345                            (double*) centroid_values -> data,
    346                            (int*) number_of_boundaries -> data,
    347                            (int*) surrogate_neighbours -> data,
    348                            (double*) a -> data,
    349                            (double*) b -> data);     
    350                            
    351   if (err != 0) {
    352     PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed");
    353     return NULL;
    354   }
    355                      
    356   //Release                 
    357   Py_DECREF(centroids);   
    358   Py_DECREF(centroid_values);   
    359   Py_DECREF(number_of_boundaries);
    360   Py_DECREF(surrogate_neighbours);           
    361          
    362   //Build result, release and return
    363   R = Py_BuildValue("OO", PyArray_Return(a), PyArray_Return(b));
    364   Py_DECREF(a);     
    365   Py_DECREF(b);       
    366   return R;
    367 }
    368 
    369 
    370 
    371 PyObject *extrapolate_second_order(PyObject *self, PyObject *args) {
    372  
    373   PyObject *quantity, *domain;
    374   PyArrayObject
    375     *centroids,            //Coordinates at centroids
    376     *centroid_values,      //Values at centroids   
    377     *vertex_coordinates,   //Coordinates at vertices           
    378     *vertex_values,        //Values at vertices       
    379     *number_of_boundaries, //Number of boundaries for each triangle
    380     *surrogate_neighbours, //True neighbours or - if one missing - self
    381     *a, *b;                //Gradients
    382    
    383   //int N, err;
    384   int dimensions[1], N, err; 
    385   //double *a, *b;  //Gradients
    386  
    387   // Convert Python arguments to C 
    388   if (!PyArg_ParseTuple(args, "O", &quantity))
    389     return NULL;
    390 
    391   domain = PyObject_GetAttrString(quantity, "domain");     
    392   if (!domain)
    393     return NULL; 
    394 
    395   //Get pertinent variables
    396   centroids = (PyArrayObject*)
    397     PyObject_GetAttrString(domain, "centroid_coordinates");
    398   if (!centroids) return NULL;   
    399  
    400   centroid_values = (PyArrayObject*)
    401     PyObject_GetAttrString(quantity, "centroid_values");   
    402   if (!centroid_values) return NULL;       
    403  
    404   surrogate_neighbours = (PyArrayObject*)
    405     PyObject_GetAttrString(domain, "surrogate_neighbours");
    406   if (!surrogate_neighbours) return NULL;       
    407  
    408   number_of_boundaries = (PyArrayObject*)
    409     PyObject_GetAttrString(domain, "number_of_boundaries");           
    410   if (!number_of_boundaries) return NULL;           
    411  
    412   vertex_coordinates = (PyArrayObject*)
    413     PyObject_GetAttrString(domain, "vertex_coordinates");           
    414   if (!vertex_coordinates) return NULL;             
    415  
    416   vertex_values = (PyArrayObject*)
    417     PyObject_GetAttrString(quantity, "vertex_values");           
    418   if (!vertex_values) return NULL;               
    419 
    420  
    421   /*
    422   printf("In extrapolate C routine\n"); 
    423   printf("d0=%d, d1=%d\n",
    424          vertex_values -> dimensions[0],
    425          vertex_values -> dimensions[1]);
    426   */
    427  
    428   vertex_values = (PyArrayObject*)
    429     PyObject_GetAttrString(quantity, "vertex_values");           
    430   if (!vertex_values) return NULL;               
    431  
    432   N = centroid_values -> dimensions[0];
    433    
    434   //Release
    435   Py_DECREF(domain);     
    436          
    437   //Allocate space for return vectors a and b (don't DECREF)
    438   dimensions[0] = N;
    439   a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
    440   b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); 
    441 
    442   //FIXME: Odd that I couldn't use normal arrays         
    443   //Allocate space for return vectors a and b (don't DECREF)
    444   //a = (double*) malloc(N * sizeof(double));
    445   //if (!a) return NULL; 
    446   //b = (double*) malloc(N * sizeof(double)); 
    447   //if (!b) return NULL;
    448 
    449  
    450   err = _compute_gradients(N,
    451                            (double*) centroids -> data,
    452                            (double*) centroid_values -> data,
    453                            (int*) number_of_boundaries -> data,
    454                            (int*) surrogate_neighbours -> data,
    455                            (double*) a -> data,
    456                            (double*) b -> data);
    457                            
    458   if (err != 0) {
    459     PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed");
    460     return NULL;
    461   }
    462  
    463   err = _extrapolate(N,
    464                      (double*) centroids -> data,
    465                      (double*) centroid_values -> data,             
    466                      (double*) vertex_coordinates -> data,
    467                      (double*) vertex_values -> data,
    468                      (double*) a -> data,
    469                      (double*) b -> data);                   
    470   //a, b);
    471                      
    472                      
    473   if (err != 0) {
    474     PyErr_SetString(PyExc_RuntimeError,
    475                     "Internal function _extrapolate failed");
    476     return NULL;
    477   }                 
    478                
    479  
    480 
    481   //Release 
    482   Py_DECREF(centroids);   
    483   Py_DECREF(centroid_values);   
    484   Py_DECREF(number_of_boundaries);
    485   Py_DECREF(surrogate_neighbours);
    486   Py_DECREF(vertex_coordinates);                   
    487   Py_DECREF(vertex_values);                 
    488   Py_DECREF(a);                   
    489   Py_DECREF(b);                     
    490   //free(a);     
    491   //free(b);       
    492  
    493   return Py_BuildValue("");   
    494 }   
    495 
    496 
    497 
    498 PyObject *limit(PyObject *self, PyObject *args) {
    499  
    500   PyObject *quantity, *domain, *Tmp;
    501   PyArrayObject
    502     *qv, //Conserved quantities at vertices
    503     *qc, //Conserved quantities at centroids
    504     *neighbours;
    505    
    506   int k, i, n, N, k3;
    507   double beta; //Safety factor
    508   double *qmin, *qmax, qn;
    509  
    510   // Convert Python arguments to C 
    511   if (!PyArg_ParseTuple(args, "O", &quantity))
    512     return NULL;
    513 
    514   domain = PyObject_GetAttrString(quantity, "domain");     
    515   if (!domain)
    516     return NULL; 
    517    
    518   neighbours = (PyArrayObject*) PyObject_GetAttrString(domain, "neighbours");
    519 
    520   //Get safety factor beta
    521   Tmp = PyObject_GetAttrString(domain, "beta");
    522   if (!Tmp)
    523     return NULL;
    524      
    525   beta = PyFloat_AsDouble(Tmp); 
    526  
    527   Py_DECREF(Tmp);   
    528   Py_DECREF(domain);     
    529 
    530   qc = (PyArrayObject*) PyObject_GetAttrString(quantity, "centroid_values");
    531   qv = (PyArrayObject*) PyObject_GetAttrString(quantity, "vertex_values");   
    532   N = qc -> dimensions[0];
    533    
    534   //Find min and max of this and neighbour's centroid values
    535   qmin = malloc(N * sizeof(double));
    536   qmax = malloc(N * sizeof(double)); 
    537   for (k=0; k<N; k++) {
    538     k3=k*3;
    539    
    540     qmin[k] = ((double*) qc -> data)[k];
    541     qmax[k] = qmin[k];
    542    
    543     for (i=0; i<3; i++) {
    544       n = ((int*) neighbours -> data)[k3+i];
    545       if (n >= 0) {
    546         qn = ((double*) qc -> data)[n]; //Neighbour's centroid value
    547          
    548         qmin[k] = min(qmin[k], qn);
    549         qmax[k] = max(qmax[k], qn);
    550       }
    551     }
    552   }
    553  
    554   // Call underlying routine
    555   _limit(N, beta, (double*) qc -> data, (double*) qv -> data, qmin, qmax);
    556          
    557   free(qmin);
    558   free(qmax); 
    559   return Py_BuildValue(""); 
    560 }
    561 
    562 
    563 
    564104// Method table for python module
    565105static struct PyMethodDef MethodTable[] = {
    566   {"limit", limit, METH_VARARGS, "Print out"},   
    567   {"update", update, METH_VARARGS, "Print out"},     
    568   {"compute_gradients", compute_gradients, METH_VARARGS, "Print out"},     
    569   {"extrapolate_second_order", extrapolate_second_order,
    570    METH_VARARGS, "Print out"},       
    571   {"interpolate_from_vertices_to_edges",
    572    interpolate_from_vertices_to_edges,
    573    METH_VARARGS, "Print out"},           
     106  {"csr_mv", csr_mv, METH_VARARGS, "Print out"},   
    574107  {NULL, NULL, 0, NULL}   /* sentinel */
    575108};
    … …  
    577110// Module initialisation   
    578111void initquantity_ext(void){
    579   Py_InitModule("quantity_ext", MethodTable);
     112  Py_InitModule("sparse_ext", MethodTable);
    580113 
    581114  import_array();     //Necessary for handling of NumPY structures 
Note: See TracChangeset for help on using the changeset viewer.