From 5da7f8935f77fa6576ba1fff35b8da42f67c5762 Mon Sep 17 00:00:00 2001 From: V-Rang Date: Wed, 3 Apr 2024 18:28:58 -0500 Subject: [PATCH] multivec added, random modified --- hippylibX/algorithms/__init__.py | 2 +- hippylibX/algorithms/multivector.py | 130 ++++++++++++++++++++++++++++ hippylibX/modeling/__init__.py | 1 + hippylibX/utils/random.py | 60 +++++++------ 4 files changed, 164 insertions(+), 29 deletions(-) create mode 100644 hippylibX/algorithms/multivector.py diff --git a/hippylibX/algorithms/__init__.py b/hippylibX/algorithms/__init__.py index 07e12147..6d96190a 100644 --- a/hippylibX/algorithms/__init__.py +++ b/hippylibX/algorithms/__init__.py @@ -1,3 +1,3 @@ from .linalg import inner # noqa from .NewtonCG import ReducedSpaceNewtonCG, ReducedSpaceNewtonCG_ParameterList # noqa -from .multivector import * #noqa \ No newline at end of file +from .multivector import * # noqa diff --git a/hippylibX/algorithms/multivector.py b/hippylibX/algorithms/multivector.py new file mode 100644 index 00000000..bae941d1 --- /dev/null +++ b/hippylibX/algorithms/multivector.py @@ -0,0 +1,130 @@ +import numpy as np +import petsc4py +import copy + + +class MultiVector: + def __init__(self, example_vec=None, nvec=None, mv=None): + if mv is None: # original + self.data = [] + self.nvec = nvec + + for i in range(self.nvec): + self.data.append(example_vec.duplicate()) + else: # copy + self.nvec = mv.nvec + self.data = copy.deepcopy(mv.data) + + def __del__(self): + for d in self.data: + d.destroy() + + def __getitem__(self, k): + return self.data[k] + + def scale(self, alpha): + if isinstance(alpha, float): + for d in self.data: + d.scale(alpha) + elif isinstance(alpha, np.array): + for i, d in enumerate(self.data): + d.scale(alpha[i]) + + def dot(self, v) -> np.array: + if isinstance(v, petsc4py.PETSc.Vec): + return_values = [] + for i in range(self.nvec): + return_values.append(self[i].dot(v)) + return np.array(return_values) + + elif isinstance(v, MultiVector): + return_values = [] + for i in range(self.nvec): + return_row = [] + for j in range(v.nvec): + return_row.append(self[i].dot(v[j])) + return_values.append(return_row) + return np.array(return_values) + + def reduce(self, alpha: np.array) -> petsc4py.PETSc.Vec: + return_vec = self[0].duplicate() + return_vec.scale(0.0) + for i in range(self.nvec): + return_vec.axpy(alpha[i], self[i]) + return return_vec + + def Borthogonalize(self, B): + return self._mgs_stable(B) + + def _mgs_stable(self, B: petsc4py.PETSc.Mat): + n = self.nvec + Bq = MultiVector(self[0], n) + r = np.zeros((n, n), dtype="d") + reorth = np.zeros((n,), dtype="d") + eps = np.finfo(np.float64).eps + + for k in np.arange(n): + B.mult(self[k], Bq[k]) + t = np.sqrt(Bq[k].dot(self[k])) + + nach = 1 + u = 0 + while nach: + u += 1 + for i in np.arange(k): + s = Bq[i].dot(self[k]) + r[i, k] += s + self[k].axpy(-s, self[i]) + + B.mult(self[k], Bq[k]) + tt = np.sqrt(Bq[k].dot(self[k])) + if tt > t * 10.0 * eps and tt < t / 10.0: + nach = 1 + t = tt + else: + nach = 0 + if tt < 10.0 * eps * t: + tt = 0.0 + + reorth[k] = u + r[k, k] = tt + if np.abs(tt * eps) > 0.0: + tt = 1.0 / tt + else: + tt = 0.0 + + self[k].scale(tt) + Bq[k].scale(tt) + + return Bq, r + + +def MatMvMult(A, x, y): + assert x.nvec() == y.nvec(), "x and y have non-matching number of vectors" + if hasattr(A, "matMvMult"): + A.matMvMult(x, y) + else: + for i in range(x.nvec()): + A.mult(x[i], y[i]) + + +def MatMvTranspmult(A, x, y): + assert x.nvec() == y.nvec(), "x and y have non-matching number of vectors" + assert hasattr(A, "transpmult"), "A does not have transpmult method implemented" + if hasattr(A, "matMvTranspmult"): + A.matMvTranspmult(x, y) + else: + for i in range(x.nvec()): + A.multTranspose(x[i], y[i]) + + +def MvDSmatMult(X, A, Y): + assert ( + X.nvec() == A.shape[0] + ), "X Number of vecs incompatible with number of rows in A" + assert ( + Y.nvec() == A.shape[1] + ), "Y Number of vecs incompatible with number of cols in A" + for j in range(Y.nvec()): + Y[j].zero() + X.reduce(Y[j], A[:, j].flatten()) diff --git a/hippylibX/modeling/__init__.py b/hippylibX/modeling/__init__.py index 6e62990e..4493c511 100644 --- a/hippylibX/modeling/__init__.py +++ b/hippylibX/modeling/__init__.py @@ -5,3 +5,4 @@ from .modelVerify import modelVerify # noqa from .Regularization import H1TikhonvFunctional, VariationalRegularization # noqa from .variables import STATE, PARAMETER, ADJOINT, NVAR # noqa +from .reducedHessian import ReducedHessian # noqa diff --git a/hippylibX/utils/random.py b/hippylibX/utils/random.py index eed94d46..b060bee9 100644 --- a/hippylibX/utils/random.py +++ b/hippylibX/utils/random.py @@ -14,43 +14,47 @@ def __init__(self, rank, nproc, seed=1): def replay(self): self.rng = np.random.MT19937(self.child_seeds[self.rank]) - def normal(self, sigma, out): - if(hasattr(out,"nvec")): - num_local_values = out[0].getLocalSize() - else: - num_local_values = out.index_map.size_local + out.index_map.num_ghosts - + def _normal_perturb_vec(self, sigma, out): + num_local_values = out.index_map.size_local + out.index_map.num_ghosts loc_random_numbers = np.random.default_rng(self.rng).normal( loc=0, scale=sigma, size=num_local_values ) - if(hasattr(out,"nvec")): #for multivec - for i in range(out.nvec): - out[i].setArray(loc_random_numbers) - else: - out.array[:] = loc_random_numbers - dlx.la.create_petsc_vector_wrap(out).ghostUpdate( - petsc4py.PETSc.InsertMode.INSERT, petsc4py.PETSc.ScatterMode.FORWARD + out.array[:] += loc_random_numbers + dlx.la.create_petsc_vector_wrap(out).ghostUpdate( + petsc4py.PETSc.InsertMode.INSERT, petsc4py.PETSc.ScatterMode.FORWARD + ) + + def _normal_perturb_multivec(self, sigma, out): + num_local_values = out[0].getLocalSize() + for i in range(out.nvec): + loc_random_numbers = np.random.default_rng(self.rng).normal( + loc=0, scale=sigma, size=num_local_values + ) + out[i].setValues( + np.arange( + out[i].getOwnershipRange()[0], + out[i].getOwnershipRange()[1], + dtype=np.int32, + ), + loc_random_numbers, + addv=petsc4py.PETSc.InsertMode.ADD, ) + out[i].assemble() def normal_perturb(self, sigma, out): - - if(hasattr(out,"nvec")): - num_local_values = out[0].getLocalSize() + if hasattr(out, "nvec"): + self._normal_perturb_multivec(sigma, out) else: - num_local_values = out.index_map.size_local + out.index_map.num_ghosts - - loc_random_numbers = np.random.default_rng(self.rng).normal( - loc=0, scale=sigma, size=num_local_values - ) - if(hasattr(out,"nvec")): #for multivec - for i in range(out.nvec): - out[i].setValues(np.arange(out[i].getOwnershipRange()[0], out[i].getOwnershipRange()[1], dtype=np.int32 ), loc_random_numbers, addv=petsc4py.PETSc.InsertMode.ADD) + self._normal_perturb_vec(sigma, out) + def normal(self, sigma, out): + if hasattr(out, "nvec"): + for d in out: + d.scale(0.0) else: - out.array[:] += loc_random_numbers - dlx.la.create_petsc_vector_wrap(out).ghostUpdate( - petsc4py.PETSc.InsertMode.INSERT, petsc4py.PETSc.ScatterMode.FORWARD - ) + out.array[:] = 0.0 + + self.normal_perturb(sigma, out) _world_rank = MPI.COMM_WORLD.rank