← Back to branch summary

~ubuntu-branches/ubuntu/maverick/dolfin/maverick

« back to all changes in this revision

Viewing changes to sandbox/la/poisson/BlockLinearAlgebra.py

  • Committer: Bazaar Package Importer
  • Author(s): Johannes Ring
  • Date: 2008-09-16 08:41:20 UTC
  • Revision ID: james.westby@ubuntu.com-20080916084120-i8k3u6lhx3mw3py3
Tags: upstream-0.9.2
Import upstream version 0.9.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
sandbox/la/poisson/BlockLinearAlgebra.py
 
1
#
 
2
# objdoc: Block Linear Algebra module
 
3
# Ola Skavhaug
 
4
# modified by Kent-Andre Mardal 
 
5
 
 
6
"""
 
7
Python block matrices and block vectors. This module provides three classes for block linear algebra systems in
 
8
Python. The main class, C{ BlockMatrix} implements a dense block matrix, 
 
9
C{DiagBlockMatrix} implements a block diagonal matrix (used e.g. for
 
10
preconditioning, and C{BlockVector} is a simple block vector consiting of
 
11
numpy.array blocks.
 
12
"""
 
13
 
 
14
import operator, numpy, math
 
15
 
 
16
class ZeroBlock(object):
 
17
    """
 
18
    This class is a zeros block matrix used to fill upp empty blocks
 
19
    in BlockMatrix.
 
20
    """
 
21
 
 
22
    def __init__(self, *args):
 
23
        """Input arguments are two integer values
 
24
        specifying the size of the block"""
 
25
        if isinstance(args[0], int) and isinstance(args[1], int):
 
26
            n = args[0]
 
27
            m = args[1]
 
28
            self.n = n
 
29
            self.m = m
 
30
        else:
 
31
            print "Wrong input value in ZeroBlock"
 
32
 
 
33
 
 
34
    def shape(self):
 
35
        """Return the number of block rows and columns."""
 
36
        return (self.n, self.m)
 
37
 
 
38
    shape = property(shape)
 
39
 
 
40
    def __mul__(self, other):
 
41
        """Matrix vector product."""
 
42
        return 0.0*other
 
43
 
 
44
 
 
45
class BlockMatrix(object):
 
46
    """A simple block matrix implementation in Python. The blocks are typically
 
47
    matrices implemented in compile languages and later exposed to Python
 
48
    either manually or using some sort of wrapper generator software."""
 
49
 
 
50
    def __init__(self, *args):
 
51
        """Input arguments are either two integer values
 
52
        specifying the number of block rows and columns, or a nested Python
 
53
        list containing the block elements directly."""
 
54
        if isinstance(args[0], int):
 
55
            n,m = args
 
56
            self.n = n
 
57
            self.m = m
 
58
            self.data = []
 
59
            self.rdim = self.cdim = None
 
60
            for i in range(n):
 
61
                self.data.append([])
 
62
                for j in range(m):
 
63
                    self.data[i].append(0)
 
64
        elif isinstance(args[0], (list, tuple)):
 
65
            self.n = len(args)
 
66
            self.m = len(args[0])
 
67
            self.data = list(args)
 
68
            self.check()
 
69
        else:
 
70
            print "ERROR", args
 
71
 
 
72
    def check(self):
 
73
        m,n = self.shape
 
74
        rdims = numpy.zeros((m,n), dtype='l')
 
75
        cdims = numpy.zeros((m,n), dtype='l')
 
76
        for i in xrange(m):
 
77
            for j in xrange(n):
 
78
                mat = self.data[i][j]
 
79
                if not isinstance(mat, int):
 
80
                    _m, _n = mat.size(0), mat.size(1)
 
81
                    rdims[i,j] = _m
 
82
                    cdims[i,j] = _n
 
83
 
 
84
        rdim = numpy.zeros(m, dtype="l")
 
85
        cdim = numpy.zeros(n, dtype="l")
 
86
        for i in xrange(m):
 
87
            row = rdims[i]
 
88
            row = numpy.where(row==0, max(row), row)
 
89
            if min(row) == max(row):
 
90
                rdim[i] = max(row)
 
91
            else:
 
92
                raise RuntimeError, "Wrong row dimensions in block matrix"
 
93
        for j in xrange(n):
 
94
            col = cdims[:,j]
 
95
            col = numpy.where(col==0, max(col), col)
 
96
            if min(col) == max(col):
 
97
                cdim[j] = max(col)
 
98
            else:
 
99
                raise RuntimeError, "Wrong col dimensions in block matrix"
 
100
        self.rdim = rdim
 
101
        self.cdim = cdim
 
102
        for i in xrange(m):
 
103
            for j in xrange(n):
 
104
                if isinstance(self.data[i][j], int):
 
105
                    import Misc
 
106
                    self.data[i][j] = Misc.genZeroBlock(self.rdim[i])
 
107
                    self.data[i][j].n = self.cdim[j]
 
108
 
 
109
            
 
110
    def row_dim(self, i):
 
111
        return self.rdim[i]
 
112
 
 
113
    def col_dim(self, j):
 
114
        return self.cdim[j]
 
115
 
 
116
    def shape(self):
 
117
        """Return the number of block rows and columns."""
 
118
        return (self.n, self.m)
 
119
 
 
120
    shape = property(shape)
 
121
 
 
122
    def __str__(self):
 
123
        """Pretty print (ugly)."""
 
124
        return str(self.data)
 
125
 
 
126
    def __repr__(self):
 
127
        """Return string capable of copying the object when calling
 
128
        C{eval(repr(self))}"""
 
129
        return "BlockMatrix(%d, %d)" % (self.n, self.m)
 
130
 
 
131
    def __setitem__(self, idx, entry):
 
132
        """Set a block matrix entry."""
 
133
        if (self._check_tuple(idx)):
 
134
            i = idx[0]; j = idx[1]
 
135
            self.data[i][j] = entry
 
136
 
 
137
    def __getitem__(self, idx):
 
138
        """Get a block matrix entry."""
 
139
        if (self._check_tuple(idx)):
 
140
            return self.data[idx[0]][idx[1]]
 
141
 
 
142
    def __add__(self, other):
 
143
        """Add two block matrices."""
 
144
        if (not isinstance(other, BlockMatrix)):
 
145
            raise TypeError, "Can not add a BlockMatrix and a %s" %(str(other))
 
146
        if (not self.shape == other.shape):
 
147
            raise ValueError, "The BlockMatrices have different shapes: %s, %s" % (str(self.shape), str(other.shape))
 
148
        b = Blockmatrix(self.n, self.m)
 
149
        for i in range(self.n):
 
150
            for j in range(self.m):
 
151
                b[i,j] = self[i]+other[j]
 
152
        return b
 
153
 
 
154
    def __mul__(self, other):
 
155
        """Matrix vector product."""
 
156
        if (not isinstance(other, BlockVector)):
 
157
            raise TypeError, "Can not multiply BlockMatrix with %s" % (str(type(other)))
 
158
        if (not other.n == self.m):
 
159
            raise ValueError, "The length of the BlockVector (%d) does not match the dimention of the BlockMatrix (%d, %d)" % (other.n, self.n, self.m)
 
160
        res = BlockVector()
 
161
        res.n = self.n
 
162
        if self.rdim == None:
 
163
            self.check()
 
164
        for i in range(self.n):
 
165
            n = self.rdim[i]
 
166
            rhs_i = other.data[i].copy()
 
167
            rhs_i.zero()
 
168
            res.data.append(rhs_i)
 
169
 
 
170
            # temporal variable
 
171
            tmp = res[i].copy()
 
172
            for j in range(self.m):
 
173
                self.data[i][j].mult(other.data[j], tmp)
 
174
                res[i].axpy(1.0,tmp)
 
175
        return res
 
176
 
 
177
    def prod(self, x, b,transposed=False): #Compute b = A*x
 
178
        """Bug in prod!!"""
 
179
        if (not isinstance(x, BlockVector)):
 
180
            raise TypeError, "Can not multiply BlockMatrix with %s" % (str(type(other)))
 
181
        if (not x.n == self.m):
 
182
            raise ValueError, "The length of the BlockVector (%d) does not match the dimention of the BlockMatrix (%d, %d)" % (x.n, self.n, self.m)
 
183
 
 
184
        for i in range(self.n):
 
185
            tmp = b[i].copy()
 
186
            tmp.zero()
 
187
            for j in range(self.m):
 
188
                if transposed: 
 
189
                    self.data[i][j].mult(x.data[j], tmp, transposed)
 
190
                else: 
 
191
                    self.data[j][i].mult(x.data[j], tmp, transposed)
 
192
#                b[i] += tmp
 
193
                b[i].axpy(1.0, tmp)
 
194
 
 
195
 
 
196
    def _check_tuple(self, idx):
 
197
        """Assure that idx is a tuple containing two-elements inside the
 
198
        range (self.n-1, self.m-1)"""
 
199
        if (not isinstance(idx, tuple)):
 
200
            raise TypeError, "wrong type: %s as index" % (str(type(idx)))
 
201
        if (not len(idx) == 2):
 
202
            raise ValueError, "length of index tuple must be 2, got %d" %(len(idx))
 
203
        if (idx[0] > self.n-1):
 
204
            raise ValueError, "index out of range, %d" %(idx[0])
 
205
        if (idx[1] > self.m-1):
 
206
            raise ValueError, "index out of range, %d" %(idx[1])
 
207
        return True
 
208
 
 
209
class DiagBlockMatrix(object):
 
210
    """A diagonal block matrix implementation in Python. This class is
 
211
    typically used to store block diagonal preconditioners"""
 
212
 
 
213
    def __init__(self, input):
 
214
        """Constructor. 
 
215
        
 
216
        Input arguments are either an integer value specifying the number of
 
217
        block rows, or a Python list containing the diagonal block elements
 
218
        directly."""
 
219
        if isinstance(input,  int):
 
220
            n = input
 
221
            self.n = n
 
222
            self.data = []
 
223
            for i in range(n):
 
224
                self.data.append([])
 
225
        elif isinstance(input, (list, tuple)):
 
226
            self.n = len(input)
 
227
            self.data = input[:]
 
228
 
 
229
 
 
230
    def shape(self):
 
231
        """Return the number of block rows."""
 
232
        return (self.n,)
 
233
 
 
234
    shape = property(shape)
 
235
 
 
236
    def __str__(self):
 
237
        """Pretty print (ugly)."""
 
238
        return str(self.data)
 
239
 
 
240
    def __repr__(self):
 
241
        """Return string capable of copying the object when calling
 
242
        eval(repr(self))"""
 
243
        return "BlockMatrix(%d)" % (self.n)
 
244
 
 
245
    def __setitem__(self, idx, entry):
 
246
        """Set a diagonal block matrix entry."""
 
247
        if (self._check_tuple(idx)):
 
248
            if not idx[0] == idx[1] :
 
249
                raise ValueError, "diagonal block matrix can not assign to indices (%d, %d)" % (idx[0], idx[1])
 
250
            self.data[idx[0]] = entry
 
251
 
 
252
    def __getitem__(self, idx):
 
253
        """Get a diagonal block matrix entry."""
 
254
        if (self._check_tuple(idx)):
 
255
            if not idx[0] == idx[1] :
 
256
                raise ValueError, "diagonal block matrix dows not provide indices (%d, %d)" % (idx[0], idx[1])
 
257
            return self.data[idx[0]]
 
258
 
 
259
    def __add__(self, other):
 
260
        """Add two diagonal block matrices."""
 
261
        if (not isinstance(other, DiagBlockMatrix)):
 
262
            raise TypeError, "Can not add a DiagBlockMatrix and a %s" %(str(other))
 
263
        if (not self.shape == other.shape):
 
264
            raise ValueError, "The BlockMatrices have different shapes: %s, %s" % (str(self.shape), str(other.shape))
 
265
        b = DiagBlockmatrix(self.n)
 
266
        for i in range(self.n):
 
267
            b.data[i] = self[i]+other[i]
 
268
        return b
 
269
 
 
270
    def __mul__(self, other):
 
271
        """Matrix vector product."""
 
272
        if (not isinstance(other, BlockVector)):
 
273
            raise TypeError, "Can not multiply DiagBlockMatrix with %s" % (str(type(other)))
 
274
        if (not other.n == self.n):
 
275
            raise ValueError, "The length of the BlockVector (%d) does not match the dimention of the DiagBlockMatrix (%d)" % (other.n, self.n)
 
276
        res = BlockVector()
 
277
        res.n = self.n
 
278
        for i in range(self.n):
 
279
            res.data.append(self.data[i]*other.data[i])
 
280
        return res
 
281
 
 
282
 
 
283
    def prod(self, x, b,transposed=False): #Compute b = A*x
 
284
        """Bug in prod!!"""
 
285
        if (not isinstance(x, BlockVector)):
 
286
            raise TypeError, "Can not multiply DiagBlockMatrix with %s" % (str(type(other)))
 
287
        if (not x.n == self.n):
 
288
            raise ValueError, "The length of the BlockVector (%d) does not match the dimention of the DiagBlockMatrix (%d, %d)" % (x.n, self.n, self.n)
 
289
        if not transposed:
 
290
            for i in range(self.n):
 
291
                b.data[i]=self.data[i]*x.data[i]
 
292
        else:
 
293
            for i in range(self.n):
 
294
                self.data[i].transposed = True
 
295
                b.data[i]=self.data[i]*x.data[i]
 
296
                self.data[i].transposed = False
 
297
#        return res
 
298
 
 
299
 
 
300
 
 
301
    def _check_tuple(self, idx):
 
302
        """Assure that idx is a tuple containing two-elements inside the
 
303
        range (self.n-1, self.m-1)"""
 
304
        if (not isinstance(idx, tuple)):
 
305
            raise TypeError, "wrong type: %s as index" % (str(type(idx)))
 
306
        if (not len(idx) == 2):
 
307
            raise ValueError, "length of index tuple must be 2, got %d" %(len(idx))
 
308
        if (idx[0] > self.n-1):
 
309
            raise ValueError, "index out of range, %d" %(idx[0])
 
310
        if (idx[1] > self.n-1):
 
311
            raise ValueError, "index out of range, %d" %(idx[1])
 
312
        return True
 
313
 
 
314
 
 
315
 
 
316
class BlockVector(object):
 
317
    """A block vector implementation in Python. 
 
318
    
 
319
    The blocks in \c BlockVector are numpy arrays. Instances of this class
 
320
    can be matrix vector multiplied with both BlockVectors and
 
321
    DiagBlockVectors."""
 
322
 
 
323
    def __init__(self, *x):
 
324
        if (x):
 
325
            self.data = list(x)
 
326
            self.n = len(x)
 
327
        else:
 
328
            self.data = []
 
329
            self.n = 0
 
330
 
 
331
    def shape(self):
 
332
        """Return the number of blocks."""
 
333
        return self.n
 
334
 
 
335
    shape = property(shape)
 
336
 
 
337
    def num_entries(self):  
 
338
        len = 0 
 
339
        for i in range(self.n):
 
340
            len += self.data[i].size()
 
341
        return len
 
342
    
 
343
 
 
344
    def norm(self):
 
345
        return reduce(operator.add, map(lambda x: numpy.dot(x,x), self.data))
 
346
 
 
347
    def inner(self, other):
 
348
        return reduce(operator.add, map(lambda x,y: x.inner(y), self.data, other.data))
 
349
 
 
350
    def __add__(self, other):
 
351
        res = BlockVector()
 
352
        res.n = self.n
 
353
        for i in range(self.n):
 
354
#            tmp = self.data[i].copy()
 
355
#            tmp.axpy(1.0, other.data[i])
 
356
#            res.data.append(tmp)
 
357
            a = self.data[i].copy()
 
358
            a += other.data[i]
 
359
            res.data.append(a)
 
360
        return res
 
361
 
 
362
    def __sub__(self, other):
 
363
        res = BlockVector()
 
364
        res.n = self.n
 
365
        for i in range(self.n):
 
366
#            tmp = self.data[i].copy()
 
367
#            tmp.axpy(-1.0, other.data[i])
 
368
#            res.data.append(tmp)
 
369
            a = self.data[i].copy()
 
370
            a -= other.data[i]
 
371
            res.data.append(a)
 
372
        return res
 
373
 
 
374
 
 
375
    def __mul__(self, other):
 
376
        res = BlockVector()
 
377
        res.n = self.n
 
378
        if isinstance(other,BlockVector):
 
379
            for i in range(self.n):
 
380
                res.data.append(self.data[i]*other[i])
 
381
        else:
 
382
            for i in range(self.n):
 
383
                tmp = self.data[i].copy()
 
384
                tmp *= other
 
385
                res.data.append(tmp)
 
386
        return res
 
387
 
 
388
    def __rmul__(self, other):
 
389
        return self.__mul__(other)
 
390
 
 
391
 
 
392
    def __imul__(self, other):
 
393
        for i in range(self.n):
 
394
            self.data[i] *= other 
 
395
        return self
 
396
 
 
397
    def __iadd__(self, other):
 
398
        for i in range(self.n):
 
399
            self.data[i] += other.data[i]
 
400
        return self
 
401
 
 
402
    def __isub__(self, other):
 
403
        for i in range(self.n):
 
404
            self.data[i] -= other.data[i]
 
405
        return self
 
406
 
 
407
    def __getitem__(self, i):
 
408
        return self.data[i]
 
409
 
 
410
    def __setitem__(self, i, x):
 
411
        if self.n > i:
 
412
            self.data[i] = x
 
413
        elif self.n == i:
 
414
            self.data.append(x)
 
415
            self.n += 1
 
416
        else:
 
417
            raise "BlockVector error","Unable to attach vector in block %d" % (i)
 
418
 
 
419
    def copy(self):
 
420
        res = BlockVector()
 
421
        res.n = self.n
 
422
        for i in range(self.n):
 
423
            res.data.append(self.data[i].copy())
 
424
        return res
 
425
 
 
426
if __name__ == '__main__':
 
427
    """Test suite"""
 
428
 
 
429
    A = BlockMatrix(2,2)
 
430
    u = numpy.arange(10,dtype='d')
 
431
    v = numpy.arange(10,dtype='d')
 
432
    x = BlockVector(u,v)
 
433
    print x.norm()
 
434
    print (x+x).norm()
 
435
    print (2*x).norm()
 
436
    x *=2
 
437
    print type(x)
 
438
    print x.norm()
 
439
    print (x+2*x).norm()
 
440
 
 
441
    B = DiagBlockMatrix(2)