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
|
#! /usr/bin/env python
import openturns as ot
import math as m
class TestHMatrixRealAssemblyFunction:
def __init__(self, vertices, scaling=1.0):
self.vertices = vertices
self.scaling = scaling
def __call__(self, i, j):
pt1 = self.vertices[i]
pt2 = self.vertices[j]
difference = pt1 - pt2
val = m.exp(-difference.norm() / self.scaling)
return val
ot.ResourceMap.SetAsBool("HMatrix-ForceSequential", True)
ot.ResourceMap.SetAsUnsignedInteger("HMatrix-MaxLeafSize", 10)
ot.PlatformInfo.SetNumericalPrecision(3)
n = 2
indices = [n, n]
intervalMesher = ot.IntervalMesher(indices)
interval = ot.Interval([0.0] * 2, [1.0] * 2)
mesh2D = intervalMesher.build(interval)
vertices = mesh2D.getVertices()
factory = ot.HMatrixFactory()
parameters = ot.HMatrixParameters()
parameters.setAssemblyEpsilon(1.0e-6)
parameters.setRecompressionEpsilon(1.0e-6)
# HMatrix must be symmetric in order to perform Cholesky decomposition
hmat = factory.build(vertices, 1, True, parameters)
simpleAssembly = TestHMatrixRealAssemblyFunction(vertices, 0.1)
hmat.assembleReal(simpleAssembly, "L")
hmatRef = ot.HMatrix(hmat)
hmat.factorize("LLt")
# Compute A - L*L^T
hmatRef.gemm("N", "T", -1.0, hmat, hmat, 1.0)
# Check LU factorization
hmat = factory.build(vertices, 1, False, parameters)
hmat.assembleReal(simpleAssembly, "N")
hmat.factorize("LU")
print("rows=", hmat.getNbRows())
print("columns=", hmat.getNbColumns())
print("norm=", ot.Point(1, hmat.norm()))
if hmatRef.norm() < 1.0e-3:
print("norm(A-LLt) < 1e-3")
else:
print("norm(A-LLt) =", hmatRef.norm())
print("diagonal=", hmat.getDiagonal())
print("compressionRatio= (%d, %d)" % hmat.compressionRatio())
print("fullrkRatio= (%d, %d)" % hmat.fullrkRatio())
# vector multiply
y = ot.Point(hmat.getNbColumns())
x = [2.0] * hmat.getNbColumns()
hmat.gemv("N", 1.0, x, 3.0, y)
print("y=", y)
# block assembly
class TestHMatrixTensorRealAssemblyFunction:
def __init__(self, covarianceModel, vertices):
self.covarianceModel = covarianceModel
self.vertices = vertices
def __call__(self, i, j):
pt1 = self.vertices[i]
pt2 = self.vertices[j]
val = self.covarianceModel(pt1, pt2)
return val
covarianceModel = ot.ExponentialModel([0.1] * 2, [1.0] * 2)
hmat = factory.build(vertices, covarianceModel.getOutputDimension(), True, parameters)
blockAssembly = TestHMatrixTensorRealAssemblyFunction(covarianceModel, vertices)
hmat.assembleTensor(blockAssembly, covarianceModel.getOutputDimension(), "L")
hmatRef = ot.HMatrix(hmat)
hmat.factorize("LLt")
normL = hmat.norm()
hmatRef.gemm("N", "T", -1.0, hmat, hmat, 1.0)
if hmatRef.norm() < 1e-3:
print("norm(A-LLt) < 1e-3")
else:
print("norm(A-LLt) =", hmatRef.norm())
alpha = 0.1
hmat.scale(alpha)
normScaled = hmat.norm()
if abs(normL - normScaled / alpha) < 1e-10:
print("|norm(L) - 10 * norm(0.1*L)| < 1e-10")
else:
print("|norm(L) - 10 * norm(0.1*L)| > 1e-10")
|
