Changeset 605
 Timestamp:
 Nov 21, 2004, 12:28:38 PM (19 years ago)
 Location:
 inundation/ga/storm_surge/pyvolution
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

inundation/ga/storm_surge/pyvolution/sparse.py
r594 r605 306 306 307 307 308 """309 Setup for C extensions310 """311 312 313 308 def csr_mv(self, B): 309 """Python version of sparse (CSR) matrix multiplication 310 """ 314 311 315 312 from Numeric import zeros, Float … … 353 350 354 351 352 353 #Setup for C extensions 355 354 import compile 356 355 if compile.can_use_C_extension('sparse_ext.c'): 357 #Replace python version with c implementations 358 from sparse_ext import csr_mv as csr_mv 359 356 #Replace python version with c implementation 357 from sparse_ext import csr_mv 360 358 361 359 if __name__ == '__main__': 
inundation/ga/storm_surge/pyvolution/sparse_ext.c
r599 r605 15 15 #include "stdio.h" 16 16 17 //Matrixvector routine 17 18 int _csr_mv(int M, 18 19 double* data, … … 24 25 long i, j, ckey; 25 26 26 //printf("test");27 27 for (i=0; i<M; i++ ) 28 28 for (ckey=row_ptr[i]; ckey<row_ptr[i+1]; ckey++) { … … 31 31 } 32 32 33 //printf("end");34 33 return 0; 35 34 } 36 35 36 //Matrixmatrix routine 37 int _csr_mm(int M, 38 int columns, 39 double* data, 40 long* colind, 41 long* row_ptr, 42 double* x, 43 double* y) { 44 45 long i, j, ckey, c, rowind_i, rowind_j; 46 47 for (i=0; i<M; i++ ) { 48 rowind_i = i*columns; 49 50 for (ckey=row_ptr[i]; ckey<row_ptr[i+1]; ckey++) { 51 j = colind[ckey]; 52 rowind_j = j*columns; 53 54 for (c=0; c<columns; c++) { 55 y[rowind_i+c] += data[ckey]*x[rowind_j+c]; 56 } 57 } 58 } 59 60 return 0; 61 } 37 62 38 63 … … 43 68 44 69 PyObject *csr_sparse, // input sparse matrix 45 *xin, *R; 70 *xin, *R; // output wrapped vector 46 71 47 72 PyArrayObject … … 53 78 54 79 55 int dimensions[1], M, err ;80 int dimensions[1], M, err, columns, rows; 56 81 57 82 // Convert Python arguments to C … … 73 98 /* printf("x.data[5] = %g\n",((double*) x>data)[5]); */ 74 99 75 100 76 101 77 102 … … 80 105 if (!data) 81 106 return NULL; 82 83 107 84 108 colind = (PyArrayObject*) … … 91 115 92 116 M = (row_ptr > dimensions[0])1; 117 118 if (x > nd == 1) { 119 // Multiplicant is a vector 120 121 //Allocate space for return vectors y (don't DECREF) 122 dimensions[0] = M; 123 y = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); 124 125 err = _csr_mv(M, 126 (double*) data > data, 127 (long*) colind > data, 128 (long*) row_ptr > data, 129 (double*) x > data, 130 (double*) y > data); 93 131 132 133 if (err != 0) { 134 PyErr_SetString(PyExc_RuntimeError, "matrix vector mult could not be calculated"); 135 return NULL; 136 } 137 } else if(x > nd == 2) { 138 94 139 95 //printf(" x.dim[0] = %i x.dim[1] = %i \n", x > dimensions[0], x > dimensions[1]); 140 rows = x > dimensions[0]; //Number of rows in x 141 columns = x > dimensions[1]; //Number of columns in x 96 142 97 //Allocate space for return vectors y (don't DECREF) 98 dimensions[0] = M; 99 y = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); 143 //Allocate space for return matrix y (don't DECREF) 144 dimensions[0] = M; //Number of rows in sparse matrix 145 dimensions[1] = columns; 146 y = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE); 147 148 err = _csr_mm(M, columns, 149 (double*) data > data, 150 (long*) colind > data, 151 (long*) row_ptr > data, 152 (double*) x > data, 153 (double*) y > data); 154 155 } else { 156 PyErr_SetString(PyExc_RuntimeError, 157 "Allowed dimensions in sparse_ext.c restricted to 1 or 2"); 158 return NULL; 159 } 100 160 101 /* for (i=0;i<25;i++){ */102 /* printf(" data[%i]=%i\n",i,(data>data)[i]); */103 /* } */104 105 /* for (i=0;i<25;i++){ */106 /* printf(" colind[%i]=%i\n",i,(colind>data)[i]); */107 /* } */108 109 err = _csr_mv(M,110 (double*) data > data,111 (long*) colind > data,112 (long*) row_ptr > data,113 (double*) x > data,114 (double*) y > data);115 116 /* for (i=0;i<25;i++){ */117 /* printf(" data[%i]=%i\n",i,data>data[i]); */118 /* } */119 120 if (err != 0) {121 PyErr_SetString(PyExc_RuntimeError, "matrix vector mult could not be calculated");122 return NULL;123 }124 161 125 162 //Release
Note: See TracChangeset
for help on using the changeset viewer.