Changeset 586
- Timestamp:
- Nov 17, 2004, 5:48:23 PM (20 years ago)
- Location:
- inundation/ga/storm_surge/pyvolution
- Files:
-
- 1 added
- 3 edited
-
sparse.py (modified) (13 diffs)
-
sparse_ext.c (added)
-
test_cg_solve.py (modified) (2 diffs)
-
test_sparse.py (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
inundation/ga/storm_surge/pyvolution/sparse.py
r485 r586 17 17 """ 18 18 19 self. A= {}19 self.Data = {} 20 20 21 21 if len(args) == 1: … … 32 32 for j in range(self.N): 33 33 if A[i, j] != 0.0: 34 self. A[i, j] = A[i, j]34 self.Data[i, j] = A[i, j] 35 35 36 36 … … 45 45 46 46 def __repr__(self): 47 return '%d X %d sparse matrix:\n' %(self.M, self.N) + `self. A`47 return '%d X %d sparse matrix:\n' %(self.M, self.N) + `self.Data` 48 48 49 49 def __len__(self): 50 50 """Return number of nonzeros of A 51 51 """ 52 return len(self. A)52 return len(self.Data) 53 53 54 54 def nonzeros(self): … … 64 64 65 65 if x != 0: 66 self. A[key] = float(x)67 else: 68 if self. A.has_key( key ):69 del self. A[key]66 self.Data[key] = float(x) 67 else: 68 if self.Data.has_key( key ): 69 del self.Data[key] 70 70 71 71 def __getitem__(self, key): … … 75 75 assert 0 <= j < self.N 76 76 77 if self. A.has_key( key ):78 return self. A[ key ]77 if self.Data.has_key( key ): 78 return self.Data[ key ] 79 79 else: 80 80 return 0.0 … … 84 84 new = Sparse(self.M,self.N) 85 85 86 for key in self. A.keys():86 for key in self.Data.keys(): 87 87 i, j = key 88 88 89 new[i,j] = self. A[i,j]89 new[i,j] = self.Data[i,j] 90 90 91 91 return new … … 99 99 for i in range(self.M): 100 100 for j in range(self.N): 101 if self. A.has_key( (i,j) ):102 D[i, j] = self. A[ (i,j) ]101 if self.Data.has_key( (i,j) ): 102 D[i, j] = self.Data[ (i,j) ] 103 103 return D 104 104 105 105 106 106 107 def __mul__(self, other): 107 108 """Multiply this matrix onto 'other' which can either be … … 129 130 130 131 #Multiply nonzero elements 131 for key in self. A.keys():132 for key in self.Data.keys(): 132 133 i, j = key 133 134 134 R[i] += self. A[key]*B[j]135 R[i] += self.Data[key]*B[j] 135 136 elif len(B.shape) == 2: 136 137 … … 142 143 #For each column 143 144 144 for key in self. A.keys():145 for key in self.Data.keys(): 145 146 i, j = key 146 147 147 R[i, col] += self. A[key]*B[j, col]148 R[i, col] += self.Data[key]*B[j, col] 148 149 149 150 … … 161 162 162 163 new = other.copy() 163 for key in self. A.keys():164 for key in self.Data.keys(): 164 165 i, j = key 165 166 166 new[i,j] += self. A[key]167 new[i,j] += self.Data[key] 167 168 168 169 return new … … 183 184 new = self.copy() 184 185 #Multiply nonzero elements 185 for key in new. A.keys():186 for key in new.Data.keys(): 186 187 i, j = key 187 188 188 new. A[key] = other*new.A[key]189 new.Data[key] = other*new.Data[key] 189 190 190 191 return new … … 213 214 214 215 #Multiply nonzero elements 215 for key in self. A.keys():216 for key in self.Data.keys(): 216 217 i, j = key 217 218 218 R[j] += self. A[key]*B[i]219 R[j] += self.Data[key]*B[i] 219 220 220 221 else: … … 223 224 return R 224 225 225 226 227 226 class Sparse_CSR: 227 228 def __init__(self, A): 229 """Create sparse matrix in csr format. 230 231 Sparse_CSR(A) #creates csr sparse matrix from sparse matrix 232 """ 233 234 from Numeric import array, Float, Int 235 236 if isinstance(A,Sparse): 237 238 from Numeric import zeros 239 keys = A.Data.keys() 240 keys.sort() 241 nnz = len(keys) 242 data = zeros ( (nnz,), Float) 243 colind = zeros ( (nnz,), Int) 244 row_ptr = zeros ( (A.M+1,), Int) 245 current_row = -1 246 k = 0 247 for key in keys: 248 ikey0 = int(key[0]) 249 ikey1 = int(key[1]) 250 if ikey0 != current_row: 251 current_row = ikey0 252 row_ptr[ikey0] = k 253 data[k] = A.Data[key] 254 colind[k] = ikey1 255 k += 1 256 for row in range(current_row+1, A.M+1): 257 row_ptr[row] = nnz 258 #row_ptr[-1] = nnz 259 260 self.data = data 261 self.colind = colind 262 self.row_ptr = row_ptr 263 self.M = A.M 264 self.N = A.N 265 else: 266 raise ValueError, "Sparse_CSR(A) expects A == Sparse Matrix" 267 268 def __repr__(self): 269 return '%d X %d sparse matrix:\n' %(self.M, self.N) + `self.data` 270 271 def __len__(self): 272 """Return number of nonzeros of A 273 """ 274 return self.row_ptr[-1] 275 276 def nonzeros(self): 277 """Return number of nonzeros of A 278 """ 279 return len(self) 280 281 def todense(self): 282 from Numeric import zeros, Float 283 284 D = zeros( (self.M, self.N), Float) 285 286 for i in range(self.M): 287 for ckey in range(self.row_ptr[i],self.row_ptr[i+1]): 288 j = self.colind[ckey] 289 D[i, j] = self.data[ckey] 290 return D 291 292 def __mul__(self, other): 293 """Multiply this matrix onto 'other' which can either be 294 a Numeric vector, a Numeric matrix or another sparse matrix. 295 """ 296 297 from Numeric import array, zeros, Float 298 299 try: 300 B = array(other) 301 except: 302 print 'FIXME: Only Numeric types implemented so far' 303 304 #Assume numeric types from now on 305 306 if len(B.shape) == 0: 307 #Scalar - use __rmul__ method 308 R = B*self 309 310 elif len(B.shape) == 1: 311 #Vector 312 assert B.shape[0] == self.N, 'Mismatching dimensions' 313 314 R = zeros(self.M, Float) #Result 315 316 #Multiply nonzero elements 317 for i in range(self.M): 318 for ckey in range(self.row_ptr[i],self.row_ptr[i+1]): 319 j = self.colind[ckey] 320 R[i] += self.data[ckey]*B[j] 321 322 elif len(B.shape) == 2: 323 324 R = zeros((self.M, B.shape[1]), Float) #Result matrix 325 326 #Multiply nonzero elements 327 328 for col in range(R.shape[1]): 329 #For each column 330 for i in range(self.M): 331 for ckey in range(self.row_ptr[i],self.row_ptr[i+1]): 332 j = self.colind[ckey] 333 R[i, col] += self.data[ckey]*B[j,col] 334 335 else: 336 raise ValueError, 'Dimension too high: d=%d' %len(B.shape) 337 338 return R 339 340 228 341 if __name__ == '__main__': 229 342 -
inundation/ga/storm_surge/pyvolution/test_cg_solve.py
r475 r586 6 6 from Numeric import dot, allclose, array, transpose, arange, ones, Float 7 7 from cg_solve import * 8 from sparse import Sparse 8 from sparse import Sparse, Sparse_CSR 9 9 10 10 … … 78 78 assert allclose(x,xe) 79 79 80 def test_solve_large_2d_csr_matrix(self): 81 """Standard 2d laplacian""" 82 83 n = 50 84 m = 50 85 86 A = Sparse(m*n, m*n) 87 88 for i in arange(0,n): 89 for j in arange(0,m): 90 I = j+m*i 91 A[I,I] = 4.0 92 if i > 0 : 93 A[I,I-m] = -1.0 94 if i < n-1 : 95 A[I,I+m] = -1.0 96 if j > 0 : 97 A[I,I-1] = -1.0 98 if j < m-1 : 99 A[I,I+1] = -1.0 100 101 xe = ones( (n*m,), Float) 102 103 # Convert to csr format 104 #print 'start covert' 105 A = Sparse_CSR(A) 106 #print 'finish covert' 107 b = A*xe 108 x = conjugate_gradient(A,b,b,iprint=10) 109 110 assert allclose(x,xe) 111 80 112 81 113 def test_solve_large_2d_with_default_guess(self): -
inundation/ga/storm_surge/pyvolution/test_sparse.py
r485 r586 164 164 C = A+B 165 165 166 assert allclose(C.todense(), [[12,0,0], [2,8,8], [0,0,4]]) 166 assert allclose(C.todense(), [[12,0,0], [2,8,8], [0,0,4]]) 167 168 def test_sparse_tocsr(self): 169 170 A = Sparse(4,3) 171 172 A[0,0] = 3 173 A[1,1] = 2 174 A[1,2] = 2 175 A[2,2] = 1 176 A[0,2] = 4 177 A[2,0] = 5 178 179 print ' ' 180 print A.todense() 181 182 B = Sparse_CSR(A) 183 184 print B.todense() 185 186 C = [1, 2, 3] 187 188 assert allclose(B*C, [15.0, 10.0 ,8.0, 0.0]) 189 190 C2 = [[1,2],[2,4],[3,6]] 191 192 print B*C2 193 194 assert allclose(B*C2, [[15.0, 30.0],[10.0, 20.0],[8.0, 16.0],[0.0, 0.0]]) 167 195 168 196
Note: See TracChangeset
for help on using the changeset viewer.
