libmatrix/laplace at master · nutils/libmatrix · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
#! /usr/bin/python

import libmatrix, traceback, numpy


def test( nprocs, size ):

  print 'testing on %d processors' % nprocs
  print 'total dofs:', size

  print 'row dof assignments:'
  bounds = ( (size-1) * numpy.arange( nprocs+1, dtype=float ) / nprocs ).astype( int )
  rowdofmap = map( numpy.arange, bounds[:-1], bounds[1:]+1 )
  print rowdofmap
  for iproc, dofs in enumerate( rowdofmap ):
    print '  proc %d:' % iproc, dofs

  comm = libmatrix.LibMatrix( nprocs=nprocs )
  comm.set_verbosity( 2 )

  params = libmatrix.ParameterList( comm )
  params['foo'] = 1
  params['bar'] = 2
  params.cprint()

  rowmap = libmatrix.Map( comm, rowdofmap, size ) # overlapping
  print 'used:', rowmap.local2global

  block = numpy.array([[1,-1],[-1,1]])
  print 'block:'
  print block

  v_npy = numpy.arange(1,size+1) * 10

  v = libmatrix.VectorBuilder( rowmap )

  for i in range(size-1):
    dofs = numpy.arange(i,i+2)
    values = v_npy[i:i+2] / 2.
    print 'adding', values, ' to ', dofs
    v.add_global( idx=[dofs], data=values )

  v = v.complete()

  v_npy[0] /= 2
  v_npy[-1] /= 2

  print v.toarray()

  print 'resulting vector:', v.toarray()
  print 'reference vector:', v_npy

  print 'resulting sum v =', v.sum()
  print 'reference sum v =', v_npy.sum()

  norm = v.norm()

  print 'resulting norm:', norm
  print 'reference norm:', numpy.linalg.norm( v_npy )

  A = libmatrix.MatrixBuilder( rowmap )
  for i in range(size-1):
    idx = numpy.arange(i,i+2)
    print 'adding block to ', idx
    A.add_global( idx=[idx,idx], data=block )
  A = A.complete()

  A_npy = numpy.eye(size,size) * 2 - numpy.eye(size,size,1) - numpy.eye(size,size,-1)
  A_npy[0,0] = A_npy[-1,-1] = 1

  print 'libmatrix norm:', A.norm()

  print 'multiplying matrix and vector'
  w = A.apply( v )

  print 'libmatrix result:', w.toarray()
  print 'reference result:', numpy.dot( A_npy, v_npy )

  B = A.toarray()
  print B
  print 'reconstructed norm:', numpy.sqrt( ( B**2 ).sum() )
  print 'reference norm:', numpy.sqrt( ( A_npy**2 ).sum() )

  rhs = libmatrix.VectorBuilder( rowmap )
  rhs.add_global( [[0]], [5] )
  rhs.add_global( [[-1]], [6] )
  rhs = rhs.complete()
  AA = libmatrix.MatrixBuilder( rowmap )
  AA.add_global( idx=[[0],[0]], data=[[1]] )
  AA.add_global( idx=[[-1],[-1]], data=[[1]] )
  AA = AA.complete()

  cons = AA.solve( rhs, symmetric=True, tol=1e-10 )
  cons.nan_from_supp( AA )

  Acons = A.constrained( cons )

  print 'Constrained matrix ='
  print Acons.toarray()
  print 'libmatrix norm:', Acons.norm()
  print 'reconstructed norm:', numpy.sqrt( ( Acons.toarray()**2 ).sum() )

  x = A.solve( precon='ILUT', constrain=cons, symmetric=True, tol=1e-10 )

  result = x.toarray()
  print 'solve result:', result
  print 'laplacian:', numpy.linalg.norm( numpy.diff( result, n=2 ) )

try:
  test( nprocs=2, size=10 )
except:
  traceback.print_exc()