From fe29b5176a2a4ae964d72b5196c31d15ec8fbc6a Mon Sep 17 00:00:00 2001 From: mateomarin97 Date: Wed, 6 Dec 2023 16:48:33 +0000 Subject: [PATCH 01/18] To array and to sparse are ready --- psydac/linalg/basic.py | 465 +++++++++++++++++- psydac/linalg/tests/test_toarray.py | 203 ++++++++ psydac/linalg/tests/test_toarray_parallel.py | 198 ++++++++ psydac/linalg/tests/test_tosparse_parallel.py | 202 ++++++++ 4 files changed, 1048 insertions(+), 20 deletions(-) create mode 100644 psydac/linalg/tests/test_toarray.py create mode 100644 psydac/linalg/tests/test_toarray_parallel.py create mode 100644 psydac/linalg/tests/test_tosparse_parallel.py diff --git a/psydac/linalg/basic.py b/psydac/linalg/basic.py index dc86957b4..29a628a75 100644 --- a/psydac/linalg/basic.py +++ b/psydac/linalg/basic.py @@ -6,6 +6,8 @@ from abc import ABC, abstractmethod from scipy.sparse import coo_matrix import numpy as np +import itertools +from scipy import sparse __all__ = ('VectorSpace', 'Vector', 'LinearOperator', 'ZeroOperator', 'IdentityOperator', 'ScaledLinearOperator', 'SumLinearOperator', 'ComposedLinearOperator', 'PowerLinearOperator', 'InverseLinearOperator', 'LinearSolver') @@ -221,14 +223,452 @@ def codomain(self): def dtype(self): pass - @abstractmethod + # Function that returns the matrix corresponding to the linear operator. Returns a scipy.sparse.csr.csr_matrix. def tosparse(self): - pass + """ + Transforms the linear operator into a matrix, which is either stored in dense or sparse format. - @abstractmethod - def toarray(self): - pass + Parameters + ---------- + out : Numpy.ndarray, optional + If given, the output will be written in-place into this array. + is_sparse : bool, optional + If set to True the method returns the matrix as a Scipy sparse matrix, if set to false + it returns the full matrix as a Numpy.ndarray + format : string, optional + Only relevant if is_sparse is True. Specifies the format in which the sparse matrix is to be stored. + Choose from "csr" (Compressed Sparse Row, default),"csc" (Compressed Sparse Column), "bsr" (Block Sparse Row ), + "lil" (List of Lists), "dok" (Dictionary of Keys), "coo" (COOrdinate format) and "dia" (DIAgonal). + Returns + ------- + out : Numpy.ndarray or scipy.sparse.csr.csr_matrix + The matrix form of the linear operator. If ran in parallel each rank gets the full + matrix representation of the linear operator. + """ + # v will be the unit vector with which we compute Av = ith column of A. + v = self.domain.zeros() + # We define a temporal vector + tmp2 = self.codomain.zeros() + + #We need to determine if we are a blockvector or a stencilvector but we are not able to use + #the BlockVectorSpace and StencilVectorSpace classes in here. So we use the following + #char. + #This char is either a b for Blockvectors or an s for Stencilvectors + BoS= str(self.domain)[15] + + if BoS == "b": + comm = self.domain.spaces[0].cart.comm + elif BoS == "s": + comm = self.domain.cart.comm + rank = comm.Get_rank() + size = comm.Get_size() + + numrows = self.codomain.dimension + numcols = self.domain.dimension + # We define a list to store the non-zero data, a list to sotre the row index of said data and a list to store the column index. + data = [] + row = [] + colarr = [] + + # V is either a BlockVector or a StencilVector depending on the domain of the linear operator. + if BoS == "b": + # we collect all starts and ends in two big lists + starts = [vi.starts for vi in v] + ends = [vi.ends for vi in v] + # We collect the dimension of the BlockVector + npts = [sp.npts for sp in self.domain.spaces] + # We get the number of space we have + nsp = len(self.domain.spaces) + # We get the number of dimensions each space has. + ndim = [sp.ndim for sp in self.domain.spaces] + + # First each rank is going to need to know the starts and ends of all other ranks + startsarr = np.array([starts[i][j] for i in range(nsp) + for j in range(ndim[i])], dtype=int) + + # Create an array to store gathered data from all ranks + allstarts = np.empty(size * len(startsarr), dtype=int) + + # Use Allgather to gather 'starts' from all ranks into 'allstarts' + comm.Allgather(startsarr, allstarts) + + # Reshape 'allstarts' to have 9 columns and 'size' rows + allstarts = allstarts.reshape((size, len(startsarr))) + + endsarr = np.array([ends[i][j] for i in range(nsp) + for j in range(ndim[i])], dtype=int) + # Create an array to store gathered data from all ranks + allends = np.empty(size * len(endsarr), dtype=int) + + # Use Allgather to gather 'ends' from all ranks into 'allends' + comm.Allgather(endsarr, allends) + + # Reshape 'allends' to have 9 columns and 'size' rows + allends = allends.reshape((size, len(endsarr))) + + currentrank = 0 + # Each rank will take care of setting to 1 each one of its entries while all other entries remain zero. + while (currentrank < size): + # since the size of npts changes denpending on h we need to compute a starting point for + # our column index + spoint = 0 + npredim = 0 + # We iterate over the stencil vectors inside the BlockVector + for h in range(nsp): + itterables = [] + for i in range(ndim[h]): + itterables.append( + range(allstarts[currentrank][i+npredim], allends[currentrank][i+npredim]+1)) + # We iterate over all the entries that belong to rank number currentrank + for i in itertools.product(*itterables): + if (rank == currentrank): + v[h][i] = 1.0 + v[h].update_ghost_regions() + # Compute dot product with the linear operator. + tmp2 *= 0. + self.dot(v, out=tmp2) + # Compute to which column this iteration belongs + col = spoint + col += np.ravel_multi_index(i, npts[h]) + aux = tmp2.toarray() + # We now need to now which entries on tmp2 are non-zero and store then in our data list + for l in np.where(aux != 0)[0]: + data.append(aux[l]) + colarr.append(col) + row.append(l) + if (rank == currentrank): + v[h][i] = 0.0 + v[h].update_ghost_regions() + cummulative = 1 + for i in range(ndim[h]): + cummulative *= npts[h][i] + spoint += cummulative + npredim += ndim[h] + currentrank += 1 + elif BoS == "s": + # We get the start and endpoint for each sublist in v + starts = v.starts + ends = v.ends + # We get the dimensions of the StencilVector + npts = self.domain.npts + # We get the number of space we have + nsp = 1 + # We get the number of dimensions the StencilVectorSpace has. + ndim = self.domain.ndim + + # First each rank is going to need to know the starts and ends of all other ranks + startsarr = np.array([starts[j] for j in range(ndim)], dtype=int) + # Create an array to store gathered data from all ranks + allstarts = np.empty(size * len(startsarr), dtype=int) + + # Use Allgather to gather 'starts' from all ranks into 'allstarts' + comm.Allgather(startsarr, allstarts) + + # Reshape 'allstarts' to have 3 columns and 'size' rows + allstarts = allstarts.reshape((size, len(startsarr))) + + endsarr = np.array([ends[j] for j in range(ndim)], dtype=int) + # Create an array to store gathered data from all ranks + allends = np.empty(size * len(endsarr), dtype=int) + + # Use Allgather to gather 'ends' from all ranks into 'allends' + comm.Allgather(endsarr, allends) + + # Reshape 'allends' to have 3 columns and 'size' rows + allends = allends.reshape((size, len(endsarr))) + + currentrank = 0 + # Each rank will take care of setting to 1 each one of its entries while all other entries remain zero. + while (currentrank < size): + itterables = [] + for i in range(ndim): + itterables.append( + range(allstarts[currentrank][i], allends[currentrank][i]+1)) + # We iterate over all the entries that belong to rank number currentrank + for i in itertools.product(*itterables): + if (rank == currentrank): + v[i] = 1.0 + v.update_ghost_regions() + # Compute dot product with the linear operator. + self.dot(v, out=tmp2) + # Compute to which column this iteration belongs + col = np.ravel_multi_index(i, npts) + aux = tmp2.toarray() + # We now need to now which entries on tmp2 are non-zero and store then in our data list + for l in np.where(aux != 0)[0]: + data.append(aux[l]) + colarr.append(col) + row.append(l) + if (rank == currentrank): + v[i] = 0.0 + v.update_ghost_regions() + currentrank += 1 + else: + # I cannot conceive any situation where this error should be thrown, but I put it here just in case something unexpected happens. + raise Exception( + 'Function toarray_struphy() only supports Stencil Vectors or Block Vectors.') + + return sparse.csr_matrix((data, (row, colarr)), shape=(numrows, numcols)) + + # Function that returns the matrix corresponding to the linear operator. Returns a numpy array. + def toarray(self, out=None, is_sparse=False, format="csr"): + """ + Transforms the linear operator into a matrix, which is either stored in dense or sparse format. + + Parameters + ---------- + out : Numpy.ndarray, optional + If given, the output will be written in-place into this array. + is_sparse : bool, optional + If set to True the method returns the matrix as a Scipy sparse matrix, if set to false + it returns the full matrix as a Numpy.ndarray + format : string, optional + Only relevant if is_sparse is True. Specifies the format in which the sparse matrix is to be stored. + Choose from "csr" (Compressed Sparse Row, default),"csc" (Compressed Sparse Column), "bsr" (Block Sparse Row ), + "lil" (List of Lists), "dok" (Dictionary of Keys), "coo" (COOrdinate format) and "dia" (DIAgonal). + + Returns + ------- + out : Numpy.ndarray or scipy.sparse.csr.csr_matrix + The matrix form of the linear operator. If ran in parallel each rank gets the full + matrix representation of the linear operator. + """ + # v will be the unit vector with which we compute Av = ith column of A. + v = self.domain.zeros() + # We define a temporal vector + tmp2 = self.codomain.zeros() + #We need to determine if we are a blockvector or a stencilvector but we are not able to use + #the BlockVectorSpace and StencilVectorSpace classes in here. So we use the following + #char. + #This char is either a b for Blockvectors or an s for Stencilvectors + BoS= str(self.domain)[15] + + if BoS == "b": + comm = self.domain.spaces[0].cart.comm + elif BoS == "s": + comm = self.domain.cart.comm + rank = comm.Get_rank() + size = comm.Get_size() + + if (is_sparse == False): + if out is None: + # We declare the matrix form of our linear operator + out = np.zeros( + [self.codomain.dimension, self.domain.dimension], dtype=self.dtype) + else: + assert isinstance(out, np.ndarray) + assert out.shape[0] == self.codomain.dimension + assert out.shape[1] == self.domain.dimension + else: + if out is not None: + raise Exception( + 'If is_sparse is True then out must be set to None.') + numrows = self.codomain.dimension + numcols = self.domain.dimension + # We define a list to store the non-zero data, a list to sotre the row index of said data and a list to store the column index. + data = [] + row = [] + colarr = [] + + # V is either a BlockVector or a StencilVector depending on the domain of the linear operator. + if BoS == "b": + # we collect all starts and ends in two big lists + starts = [vi.starts for vi in v] + ends = [vi.ends for vi in v] + # We collect the dimension of the BlockVector + npts = [sp.npts for sp in self.domain.spaces] + # We get the number of space we have + nsp = len(self.domain.spaces) + # We get the number of dimensions each space has. + ndim = [sp.ndim for sp in self.domain.spaces] + + # First each rank is going to need to know the starts and ends of all other ranks + startsarr = np.array([starts[i][j] for i in range(nsp) + for j in range(ndim[i])], dtype=int) + + # Create an array to store gathered data from all ranks + allstarts = np.empty(size * len(startsarr), dtype=int) + + # Use Allgather to gather 'starts' from all ranks into 'allstarts' + comm.Allgather(startsarr, allstarts) + + # Reshape 'allstarts' to have 9 columns and 'size' rows + allstarts = allstarts.reshape((size, len(startsarr))) + + endsarr = np.array([ends[i][j] for i in range(nsp) + for j in range(ndim[i])], dtype=int) + # Create an array to store gathered data from all ranks + allends = np.empty(size * len(endsarr), dtype=int) + + # Use Allgather to gather 'ends' from all ranks into 'allends' + comm.Allgather(endsarr, allends) + + # Reshape 'allends' to have 9 columns and 'size' rows + allends = allends.reshape((size, len(endsarr))) + + currentrank = 0 + # Each rank will take care of setting to 1 each one of its entries while all other entries remain zero. + while (currentrank < size): + # since the size of npts changes denpending on h we need to compute a starting point for + # our column index + spoint = 0 + npredim = 0 + # We iterate over the stencil vectors inside the BlockVector + for h in range(nsp): + itterables = [] + for i in range(ndim[h]): + itterables.append( + range(allstarts[currentrank][i+npredim], allends[currentrank][i+npredim]+1)) + # We iterate over all the entries that belong to rank number currentrank + for i in itertools.product(*itterables): + if (rank == currentrank): + v[h][i] = 1.0 + v[h].update_ghost_regions() + # Compute dot product with the linear operator. + tmp2 *= 0. + self.dot(v, out=tmp2) + # Compute to which column this iteration belongs + col = spoint + col += np.ravel_multi_index(i, npts[h]) + if is_sparse == False: + out[:, col] = tmp2.toarray() + else: + aux = tmp2.toarray() + # We now need to now which entries on tmp2 are non-zero and store then in our data list + for l in np.where(aux != 0)[0]: + data.append(aux[l]) + colarr.append(col) + row.append(l) + if (rank == currentrank): + v[h][i] = 0.0 + v[h].update_ghost_regions() + cummulative = 1 + for i in range(ndim[h]): + cummulative *= npts[h][i] + spoint += cummulative + npredim += ndim[h] + currentrank += 1 + elif BoS == "s": + # We get the start and endpoint for each sublist in v + starts = v.starts + ends = v.ends + # We get the dimensions of the StencilVector + npts = self.domain.npts + # We get the number of space we have + nsp = 1 + # We get the number of dimensions the StencilVectorSpace has. + ndim = self.domain.ndim + + # First each rank is going to need to know the starts and ends of all other ranks + startsarr = np.array([starts[j] for j in range(ndim)], dtype=int) + # Create an array to store gathered data from all ranks + allstarts = np.empty(size * len(startsarr), dtype=int) + + # Use Allgather to gather 'starts' from all ranks into 'allstarts' + comm.Allgather(startsarr, allstarts) + + # Reshape 'allstarts' to have 3 columns and 'size' rows + allstarts = allstarts.reshape((size, len(startsarr))) + + endsarr = np.array([ends[j] for j in range(ndim)], dtype=int) + # Create an array to store gathered data from all ranks + allends = np.empty(size * len(endsarr), dtype=int) + + # Use Allgather to gather 'ends' from all ranks into 'allends' + comm.Allgather(endsarr, allends) + + # Reshape 'allends' to have 3 columns and 'size' rows + allends = allends.reshape((size, len(endsarr))) + + currentrank = 0 + # Each rank will take care of setting to 1 each one of its entries while all other entries remain zero. + while (currentrank < size): + itterables = [] + for i in range(ndim): + itterables.append( + range(allstarts[currentrank][i], allends[currentrank][i]+1)) + # We iterate over all the entries that belong to rank number currentrank + for i in itertools.product(*itterables): + if (rank == currentrank): + v[i] = 1.0 + v.update_ghost_regions() + # Compute dot product with the linear operator. + self.dot(v, out=tmp2) + # Compute to which column this iteration belongs + col = np.ravel_multi_index(i, npts) + if is_sparse == False: + out[:, col] = tmp2.toarray() + else: + aux = tmp2.toarray() + # We now need to now which entries on tmp2 are non-zero and store then in our data list + for l in np.where(aux != 0)[0]: + data.append(aux[l]) + colarr.append(col) + row.append(l) + if (rank == currentrank): + v[i] = 0.0 + v.update_ghost_regions() + currentrank += 1 + else: + # I cannot conceive any situation where this error should be thrown, but I put it here just in case something unexpected happens. + raise Exception( + 'Function toarray_struphy() only supports Stencil Vectors or Block Vectors.') + + if is_sparse == False: + return out + else: + # Gather all rows on rank 0 + gathered_rows = comm.gather(row, root=0) + # Gather all colarr on rank 0 + gathered_cols = comm.gather(colarr, root=0) + # Gather all data on rank 0 + gathered_data = comm.gather(data, root=0) + + if rank == 0: + # Rank 0 collects all rows from other ranks + all_rows = [ + item for sublist in gathered_rows for item in sublist] + # Rank 0 collects all columns from other ranks + all_cols = [ + item for sublist in gathered_cols for item in sublist] + # Rank 0 collects all data from other ranks + all_data = [ + item for sublist in gathered_data for item in sublist] + + # Broadcast 'all_rows' to all other ranks + comm.bcast(all_rows, root=0) + # Broadcast 'all_cols' to all other ranks + comm.bcast(all_cols, root=0) + # Broadcast 'all_data' to all other ranks + comm.bcast(all_data, root=0) + else: + # Other ranks receive the 'all_rows' list through broadcast + all_rows = comm.bcast(None, root=0) + # Other ranks receive the 'all_cols' list through broadcast + all_cols = comm.bcast(None, root=0) + # Other ranks receive the 'all_data' list through broadcast + all_data = comm.bcast(None, root=0) + + if format == "csr": + return sparse.csr_matrix((all_data, (all_rows, all_cols)), shape=(numrows, numcols)) + elif format == "csc": + return sparse.csc_matrix((all_data, (all_rows, all_cols)), shape=(numrows, numcols)) + elif format == "bsr": + return sparse.bsr_matrix((all_data, (all_rows, all_cols)), shape=(numrows, numcols)) + elif format == "lil": + return sparse.csr_matrix((all_data, (all_rows, all_cols)), shape=(numrows, numcols)).tolil() + elif format == "dok": + return sparse.csr_matrix((all_data, (all_rows, all_cols)), shape=(numrows, numcols)).todok() + elif format == "coo": + return sparse.coo_matrix((all_data, (all_rows, all_cols)), shape=(numrows, numcols)) + elif format == "dia": + return sparse.csr_matrix((all_data, (all_rows, all_cols)), shape=(numrows, numcols)).todia() + else: + raise Exception( + 'The selected sparse matrix format must be one of the following : csr, csc, bsr, lil, dok, coo or dia.') + + @abstractmethod def dot(self, v, out=None): """ Apply linear operator to Vector v. Result is written to Vector out, if provided.""" @@ -748,9 +1188,6 @@ def multiplicants(self): def dtype(self): return None - def toarray(self): - raise NotImplementedError('toarray() is not defined for ComposedLinearOperators.') - def tosparse(self): mats = [M.tosparse() for M in self._multiplicants] M = mats[0] @@ -849,12 +1286,6 @@ def operator(self): def factorial(self): return self._factorial - def toarray(self): - raise NotImplementedError('toarray() is not defined for PowerLinearOperators.') - - def tosparse(self): - raise NotImplementedError('tosparse() is not defined for PowerLinearOperators.') - def transpose(self, conjugate=False): return PowerLinearOperator(domain=self._codomain, codomain=self._domain, A=self._operator.transpose(conjugate=conjugate), n=self._factorial) @@ -905,12 +1336,6 @@ def linop(self): def options(self): return self._options - def toarray(self): - raise NotImplementedError('toarray() is not defined for InverseLinearOperators.') - - def tosparse(self): - raise NotImplementedError('tosparse() is not defined for InverseLinearOperators.') - def get_info(self): return self._info diff --git a/psydac/linalg/tests/test_toarray.py b/psydac/linalg/tests/test_toarray.py new file mode 100644 index 000000000..f347613ce --- /dev/null +++ b/psydac/linalg/tests/test_toarray.py @@ -0,0 +1,203 @@ +import pytest +import numpy as np + +from psydac.linalg.block import BlockLinearOperator, BlockVector, BlockVectorSpace +from psydac.linalg.basic import LinearOperator, ZeroOperator, IdentityOperator, ComposedLinearOperator, InverseLinearOperator, SumLinearOperator, PowerLinearOperator, ScaledLinearOperator +from psydac.linalg.stencil import StencilVectorSpace, StencilVector, StencilMatrix +from psydac.linalg.solvers import ConjugateGradient, inverse +from psydac.ddm.cart import DomainDecomposition, CartDecomposition +from mpi4py import MPI +from random import random, seed +#=============================================================================== + +n1array = [2, 7] +n2array = [2, 3] +p1array = [1, 3] +p2array = [1, 3] + + +def compute_global_starts_ends(domain_decomposition, npts): + ndims = len(npts) + global_starts = [None]*ndims + global_ends = [None]*ndims + + for axis in range(ndims): + es = domain_decomposition.global_element_starts[axis] + ee = domain_decomposition.global_element_ends [axis] + + global_ends [axis] = ee.copy() + global_ends [axis][-1] = npts[axis]-1 + global_starts[axis] = np.array([0] + (global_ends[axis][:-1]+1).tolist()) + + return global_starts, global_ends + +def get_StencilVectorSpace(n1, n2, p1, p2, P1, P2): + comm = MPI.COMM_WORLD + npts = [n1, n2] + pads = [p1, p2] + periods = [P1, P2] + D = DomainDecomposition(npts, periods=periods, comm = comm) + global_starts, global_ends = compute_global_starts_ends(D, npts) + C = CartDecomposition(D, npts, global_starts, global_ends, pads=pads, shifts=[1,1]) + V = StencilVectorSpace(C) + return V + +#=============================================================================== +# SERIAL TESTS +#=============================================================================== +@pytest.mark.parametrize('n1', n1array) +@pytest.mark.parametrize('n2', n2array) +@pytest.mark.parametrize('p1', p1array) +@pytest.mark.parametrize('p2', p2array) + +def test_square_stencil_basic(n1, n2, p1, p2, P1=False, P2=False): + + # 1. Initiate square LOs S,S1 (StencilMatrix), I (IdentityOperator), Z (ZeroOperator) and a Stencilvector v + # 2. Test general basic operations + # 3. Test special cases + + ### + ### 1. Initiation + ### + + # Initiate StencilVectorSpace + V = get_StencilVectorSpace(n1, n2, p1, p2, P1, P2) + + S = StencilMatrix(V, V) + + # Initiate a StencilVector + v = StencilVector(V) + for i in range(n1): + for j in range(n2): + v[i,j] = 1 + + nonzero_values = dict() + for k1 in range(-p1,p1+1): + for k2 in range(-p2,p2+1): + nonzero_values[k1,k2] = 1 + k1*n2 + k2 + for k1 in range(-p1,p1+1): + for k2 in range(-p2,p2+1): + if k1==0: + if k2<0: + nonzero_values[k1,k2] = nonzero_values[-k1,-k2] + elif k1<0: + nonzero_values[k1,k2] = nonzero_values[-k1,-k2] + + for k1 in range(-p1,p1+1): + for k2 in range(-p2,p2+1): + S[:,:,k1,k2] = nonzero_values[k1,k2] + S.remove_spurious_entries() + + # Construct exact matrices by hand + A1 = np.zeros( S.shape ) + for i1 in range(n1): + for i2 in range(n2): + for k1 in range(-p1,p1+1): + for k2 in range(-p2,p2+1): + j1 = (i1+k1) % n1 + j2 = (i2+k2) % n2 + i = i1*(n2) + i2 + j = j1*(n2) + j2 + if (P1 or 0 <= i1+k1 < n1) and (P2 or 0 <= i2+k2 < n2): + A1[i,j] = nonzero_values[k1,k2] + + + ### + ### Use PowerLinearOperator to test toarray function. + ### + ### + # Raising a StencilMatrix to a power works + Spo3 = S**3 + + #We call toarray on the PowerLinearOperator + Spo3arr = Spo3.toarray() + + #Again we know that A1 is the matrix representation of S. Thus, if we take A1 A1 A1 it should be the same as Sp3arr + Apoarr = np.matmul(A1, A1) + Apoarr = np.matmul(Apoarr, A1) + + assert np.array_equal(Spo3arr, Apoarr) + + ### + ### Use InverseLinearOperator to test toarray function. + ### + ### + #We invert S + Sinv = inverse(S, "bicgstab") + Sinvarr = Sinv.toarray() + + assert isinstance(Sinv, InverseLinearOperator) + + Itest = np.matmul(A1,Sinvarr) + assert np.allclose(Itest, np.eye(np.size(Itest[0])),atol= 10**(-5)) + + +#=============================================================================== +@pytest.mark.parametrize('n1', n1array) +@pytest.mark.parametrize('n2', n2array) +@pytest.mark.parametrize('p1', p1array) +@pytest.mark.parametrize('p2', p2array) + +def test_square_block_basic(n1, n2, p1, p2, P1=False, P2=False): + + # 1. Initiate square LOs S,S1 (StencilMatrix), Z (ZeroOperator) and a Stencilvector v + # Initiate square LOs B,B1 (BlockLO), BZ (ZeroOperator), BI (IdentityOperator) and a BlockVector vb + # 2. Test general basic operations + # 3. Test special cases + + # Initiate StencilVectorSpace + V = get_StencilVectorSpace(n1, n2, p1, p2, P1, P2) + + # Initiate Linear Operators + S = StencilMatrix(V, V) + + # Initiate a StencilVector + v = StencilVector(V) + for i in range(n1): + for j in range(n2): + v[i,j] = float(i*3.5 + j*10.0) + + nonzero_values = dict() + for k1 in range(-p1,p1+1): + for k2 in range(-p2,p2+1): + nonzero_values[k1,k2] = 1 + k1*n2 + k2 + for k1 in range(-p1,p1+1): + for k2 in range(-p2,p2+1): + if k1==0: + if k2<0: + nonzero_values[k1,k2] = nonzero_values[-k1,-k2] + elif k1<0: + nonzero_values[k1,k2] = nonzero_values[-k1,-k2] + + for k1 in range(-p1,p1+1): + for k2 in range(-p2,p2+1): + S[:,:,k1,k2] = nonzero_values[k1,k2] + S.remove_spurious_entries() + + # Initiate a BlockVectorSpace + Vb = BlockVectorSpace(V,V) + + # Initiate BlockLOs and LOs acting on BlockVectorSpaces + B = BlockLinearOperator(Vb, Vb, ((S, None), (None, S))) + + # Initiate a BlockVector + vb = BlockVector(Vb, (v, v)) + vbarr = vb.toarray() + + ### + ### Test toarray with power of one BlockLinearOperators and taking dot product with a block vector + ### + ### + Bpow = B**3 + Bpowarr = Bpow.toarray() + assert isinstance(Bpow, PowerLinearOperator) + assert np.array_equal(np.matmul(Bpowarr,vbarr),Bpow.dot(vb).toarray()) + + + +#=============================================================================== +# SCRIPT FUNCTIONALITY +#=============================================================================== +if __name__ == "__main__": + import sys + pytest.main( sys.argv ) \ No newline at end of file diff --git a/psydac/linalg/tests/test_toarray_parallel.py b/psydac/linalg/tests/test_toarray_parallel.py new file mode 100644 index 000000000..9580f3cfe --- /dev/null +++ b/psydac/linalg/tests/test_toarray_parallel.py @@ -0,0 +1,198 @@ +# -*- coding: UTF-8 -*- +# +import pytest +import numpy as np +import scipy.sparse as spa + +from psydac.linalg.basic import LinearOperator, ZeroOperator, IdentityOperator, ComposedLinearOperator, InverseLinearOperator, SumLinearOperator, PowerLinearOperator, ScaledLinearOperator +from psydac.linalg.direct_solvers import SparseSolver +from psydac.linalg.stencil import StencilVectorSpace, StencilVector, StencilMatrix +from psydac.linalg.block import BlockVectorSpace, BlockVector +from psydac.linalg.block import BlockLinearOperator +from psydac.linalg.utilities import array_to_psydac +from psydac.linalg.kron import KroneckerLinearSolver +from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL +from psydac.ddm.cart import DomainDecomposition, CartDecomposition + +#=============================================================================== +def compare_arrays(arr_psy, arr, rank, atol=1e-14, verbose=False): + '''Assert equality of distributed psydac array and corresponding fraction of cloned numpy array. + Arrays can be block-structured as nested lists/tuples. + + Parameters + ---------- + arr_psy : psydac object + Stencil/Block Vector/Matrix. + + arr : array + Numpy array with same tuple/list structure as arr_psy. If arr is a matrix it can be stencil or band format. + + rank : int + Rank of mpi process + + atol : float + Absolute tolerance used in numpy.allclose() + ''' + if isinstance(arr_psy, StencilVector): + + s = arr_psy.starts + e = arr_psy.ends + + tmp1 = arr_psy[s[0]: e[0] + 1, s[1]: e[1] + 1] + + if arr.ndim == 2: + tmp2 = arr[s[0]: e[0] + 1, s[1]: e[1] + 1] + else: + tmp2 = arr.reshape(arr_psy.space.npts[0], + arr_psy.space.npts[1])[s[0]: e[0] + 1, s[1]: e[1] + 1] + + assert np.allclose(tmp1, tmp2, atol=atol) + + elif isinstance(arr_psy, BlockVector): + + if not (isinstance(arr, tuple) or isinstance(arr, list)): + arrs = np.split(arr, [arr_psy.blocks[0].shape[0], + arr_psy.blocks[0].shape[0] + arr_psy.blocks[1].shape[0]]) + else: + arrs = arr + + for vec_psy, vec in zip(arr_psy.blocks, arrs): + s = vec_psy.starts + e = vec_psy.ends + + tmp1 = vec_psy[s[0]: e[0] + 1, s[1]: e[1] + 1] + + if vec.ndim == 2: + tmp2 = vec[s[0]: e[0] + 1, s[1]: e[1] + 1] + else: + tmp2 = vec.reshape(vec_psy.space.npts[0], vec_psy.space.npts[1])[ + s[0]: e[0] + 1, s[1]: e[1] + 1] + + assert np.allclose(tmp1, tmp2, atol=atol) + else: + raise AssertionError('Wrong input type.') + + if verbose: + print( + f'Rank {rank}: Assertion for array comparison passed with atol={atol}.') + +def compute_global_starts_ends(domain_decomposition, npts): + global_starts = [None]*len(npts) + global_ends = [None]*len(npts) + + for axis in range(len(npts)): + es = domain_decomposition.global_element_starts[axis] + ee = domain_decomposition.global_element_ends [axis] + + global_ends [axis] = ee.copy() + global_ends [axis][-1] = npts[axis]-1 + global_starts[axis] = np.array([0] + (global_ends[axis][:-1]+1).tolist()) + + return global_starts, global_ends + +def create_equal_random_arrays(W, seedv =123): + np.random.seed(seedv) + arr = [] + arr_psy = BlockVector(W) + + for d, block in enumerate(arr_psy.blocks): + + dims = W.spaces[d].npts + + arr += [np.random.rand(*dims)] + + s = block.starts + e = block.ends + + arr_psy[d][s[0]:e[0] + 1, s[1]:e[1] + 1] = arr[-1][s[0]:e[0] + 1, s[1]:e[1] + 1] + arr_psy.update_ghost_regions() + + return arr, arr_psy + + + +@pytest.mark.parametrize( 'dtype', [float] ) +@pytest.mark.parametrize( 'n1', [8, 16] ) +@pytest.mark.parametrize( 'n2', [8, 32] ) +@pytest.mark.parametrize( 'p1', [1, 3] ) +@pytest.mark.parametrize( 'p2', [2] ) +@pytest.mark.parametrize( 'P1', [True, False] ) +@pytest.mark.parametrize( 'P2', [True] ) +@pytest.mark.parallel + +def test_block_linear_operator_parallel_dot( dtype, n1, n2, p1, p2, P1, P2 ): + # Define a factor for the data + if dtype == complex: + factor = 1j + else: + factor = 1 + + # set seed for reproducibility + np.random.seed(n1*n2*p1*p2) + + from mpi4py import MPI + + comm = MPI.COMM_WORLD + rank = comm.Get_rank() + size = comm.Get_size() + + D = DomainDecomposition([n1,n2], periods=[P1,P2], comm=comm) + + # Partition the points + npts = [n1,n2] + global_starts, global_ends = compute_global_starts_ends(D, npts) + + cart = CartDecomposition(D, npts, global_starts, global_ends, pads=[p1,p2], shifts=[1,1]) + + # Create vector space, stencil matrix, and stencil vector + V = StencilVectorSpace( cart, dtype=dtype ) + # Create and Fill Block objects + W = BlockVectorSpace(V, V) + + v1arr1, v1 = create_equal_random_arrays(W, seedv=4568) + v1arr = [] + for i in v1arr1: + aux = i.flatten() + for j in aux: + v1arr.append(j) + + # Test copy with an out + # Create random matrix + N1 = StencilMatrix( V, V ) + N2 = StencilMatrix( V, V ) + N3 = StencilMatrix( V, V ) + N4 = StencilMatrix( V, V ) + + for k1 in range(-p1,p1+1): + for k2 in range(-p2,p2+1): + N1[:,:,k1,k2] = factor*np.random.random() + N2[:,:,k1,k2] = factor*np.random.random() + N3[:,:,k1,k2] = factor*np.random.random() + N4[:,:,k1,k2] = factor*np.random.random() + + K = BlockLinearOperator( W, W ) + + K[0,0] = N1 + K[0,1] = N2 + K[1,0] = N3 + K[1,1] = N4 + + #Now we use a PowerLinearOperator + KP = K**3 + KParrloc = KP.toarray() + KParr = np.zeros(np.shape(KParrloc)) + comm.Allreduce(KParrloc, KParr, op=MPI.SUM) + + + assert isinstance(KP, PowerLinearOperator) + compare_arrays(KP.dot(v1), np.matmul(KParr, v1arr), rank) + + + + +#=============================================================================== +# SCRIPT FUNCTIONALITY +#=============================================================================== +if __name__ == "__main__": + import sys + pytest.main( sys.argv ) \ No newline at end of file diff --git a/psydac/linalg/tests/test_tosparse_parallel.py b/psydac/linalg/tests/test_tosparse_parallel.py new file mode 100644 index 000000000..4eb4a5c46 --- /dev/null +++ b/psydac/linalg/tests/test_tosparse_parallel.py @@ -0,0 +1,202 @@ +# -*- coding: UTF-8 -*- +# +import pytest +import numpy as np +import scipy.sparse as spa + +from psydac.linalg.basic import LinearOperator, ZeroOperator, IdentityOperator, ComposedLinearOperator, InverseLinearOperator, SumLinearOperator, PowerLinearOperator, ScaledLinearOperator +from psydac.linalg.direct_solvers import SparseSolver +from psydac.linalg.stencil import StencilVectorSpace, StencilVector, StencilMatrix +from psydac.linalg.block import BlockVectorSpace, BlockVector +from psydac.linalg.block import BlockLinearOperator +from psydac.linalg.utilities import array_to_psydac +from psydac.linalg.kron import KroneckerLinearSolver +from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL +from psydac.ddm.cart import DomainDecomposition, CartDecomposition + +#=============================================================================== +def compare_arrays(arr_psy, arr, rank, atol=1e-14, verbose=False): + '''Assert equality of distributed psydac array and corresponding fraction of cloned numpy array. + Arrays can be block-structured as nested lists/tuples. + + Parameters + ---------- + arr_psy : psydac object + Stencil/Block Vector/Matrix. + + arr : array + Numpy array with same tuple/list structure as arr_psy. If arr is a matrix it can be stencil or band format. + + rank : int + Rank of mpi process + + atol : float + Absolute tolerance used in numpy.allclose() + ''' + if isinstance(arr_psy, StencilVector): + + s = arr_psy.starts + e = arr_psy.ends + + tmp1 = arr_psy[s[0]: e[0] + 1, s[1]: e[1] + 1] + + if arr.ndim == 2: + tmp2 = arr[s[0]: e[0] + 1, s[1]: e[1] + 1] + else: + tmp2 = arr.reshape(arr_psy.space.npts[0], + arr_psy.space.npts[1])[s[0]: e[0] + 1, s[1]: e[1] + 1] + + assert np.allclose(tmp1, tmp2, atol=atol) + + elif isinstance(arr_psy, BlockVector): + + if not (isinstance(arr, tuple) or isinstance(arr, list)): + arrs = np.split(arr, [arr_psy.blocks[0].shape[0], + arr_psy.blocks[0].shape[0] + arr_psy.blocks[1].shape[0]]) + else: + arrs = arr + + for vec_psy, vec in zip(arr_psy.blocks, arrs): + s = vec_psy.starts + e = vec_psy.ends + + tmp1 = vec_psy[s[0]: e[0] + 1, s[1]: e[1] + 1] + + if vec.ndim == 2: + tmp2 = vec[s[0]: e[0] + 1, s[1]: e[1] + 1] + else: + tmp2 = vec.reshape(vec_psy.space.npts[0], vec_psy.space.npts[1])[ + s[0]: e[0] + 1, s[1]: e[1] + 1] + + assert np.allclose(tmp1, tmp2, atol=atol) + else: + raise AssertionError('Wrong input type.') + + if verbose: + print( + f'Rank {rank}: Assertion for array comparison passed with atol={atol}.') + +def compute_global_starts_ends(domain_decomposition, npts): + global_starts = [None]*len(npts) + global_ends = [None]*len(npts) + + for axis in range(len(npts)): + es = domain_decomposition.global_element_starts[axis] + ee = domain_decomposition.global_element_ends [axis] + + global_ends [axis] = ee.copy() + global_ends [axis][-1] = npts[axis]-1 + global_starts[axis] = np.array([0] + (global_ends[axis][:-1]+1).tolist()) + + return global_starts, global_ends + +def create_equal_random_arrays(W, seedv =123): + np.random.seed(seedv) + arr = [] + arr_psy = BlockVector(W) + + for d, block in enumerate(arr_psy.blocks): + + dims = W.spaces[d].npts + + arr += [np.random.rand(*dims)] + + s = block.starts + e = block.ends + + arr_psy[d][s[0]:e[0] + 1, s[1]:e[1] + 1] = arr[-1][s[0]:e[0] + 1, s[1]:e[1] + 1] + arr_psy.update_ghost_regions() + + return arr, arr_psy + + + +@pytest.mark.parametrize( 'dtype', [float] ) +@pytest.mark.parametrize( 'n1', [8, 16] ) +@pytest.mark.parametrize( 'n2', [8, 32] ) +@pytest.mark.parametrize( 'p1', [1, 3] ) +@pytest.mark.parametrize( 'p2', [2] ) +@pytest.mark.parametrize( 'P1', [True, False] ) +@pytest.mark.parametrize( 'P2', [True] ) +@pytest.mark.parallel + +def test_block_linear_operator_parallel_dot( dtype, n1, n2, p1, p2, P1, P2 ): + # Define a factor for the data + if dtype == complex: + factor = 1j + else: + factor = 1 + + # set seed for reproducibility + np.random.seed(n1*n2*p1*p2) + + from mpi4py import MPI + + comm = MPI.COMM_WORLD + rank = comm.Get_rank() + size = comm.Get_size() + + D = DomainDecomposition([n1,n2], periods=[P1,P2], comm=comm) + + # Partition the points + npts = [n1,n2] + global_starts, global_ends = compute_global_starts_ends(D, npts) + + cart = CartDecomposition(D, npts, global_starts, global_ends, pads=[p1,p2], shifts=[1,1]) + + # Create vector space, stencil matrix, and stencil vector + V = StencilVectorSpace( cart, dtype=dtype ) + + # Create and Fill Block objects + W = BlockVectorSpace(V, V) + + v1arr1, v1 = create_equal_random_arrays(W, seedv=4568) + v1arr = [] + for i in v1arr1: + aux = i.flatten() + for j in aux: + v1arr.append(j) + + # Test copy with an out + # Create random matrix + N1 = StencilMatrix( V, V ) + N2 = StencilMatrix( V, V ) + N3 = StencilMatrix( V, V ) + N4 = StencilMatrix( V, V ) + + for k1 in range(-p1,p1+1): + for k2 in range(-p2,p2+1): + N1[:,:,k1,k2] = factor*np.random.random() + N2[:,:,k1,k2] = factor*np.random.random() + N3[:,:,k1,k2] = factor*np.random.random() + N4[:,:,k1,k2] = factor*np.random.random() + + K = BlockLinearOperator( W, W ) + + K[0,0] = N1 + K[0,1] = N2 + K[1,0] = N3 + K[1,1] = N4 + + #We take a power of the LinearOperator + KP = K**3 + #We turn our PowerLinearOperator into a sparse matrix, and then into a numpy array + KLarrloc = KP.tosparse().toarray() + #We get the global matrix of the operator + KLarr = np.zeros(np.shape(KLarrloc)) + comm.Allreduce(KLarrloc, KLarr, op=MPI.SUM) + + assert isinstance(KP, PowerLinearOperator) + compare_arrays(KP.dot(v1), KLarr.dot(v1arr), rank) + + + + + + +#=============================================================================== +# SCRIPT FUNCTIONALITY +#=============================================================================== +if __name__ == "__main__": + import sys + pytest.main( sys.argv ) \ No newline at end of file From 8deac4aed9e14d3a52a2ba80e5f6cf2543665c22 Mon Sep 17 00:00:00 2001 From: mateomarin97 Date: Wed, 6 Dec 2023 17:02:55 +0000 Subject: [PATCH 02/18] Cleaned the comments --- psydac/linalg/basic.py | 140 +++--------------- psydac/linalg/tests/test_toarray.py | 10 +- psydac/linalg/tests/test_toarray_parallel.py | 5 +- psydac/linalg/tests/test_tosparse_parallel.py | 3 +- 4 files changed, 26 insertions(+), 132 deletions(-) diff --git a/psydac/linalg/basic.py b/psydac/linalg/basic.py index 29a628a75..e31d975a2 100644 --- a/psydac/linalg/basic.py +++ b/psydac/linalg/basic.py @@ -223,27 +223,15 @@ def codomain(self): def dtype(self): pass - # Function that returns the matrix corresponding to the linear operator. Returns a scipy.sparse.csr.csr_matrix. + # Function that returns the local matrix corresponding to the linear operator. Returns a scipy.sparse.csr.csr_matrix. def tosparse(self): """ - Transforms the linear operator into a matrix, which is either stored in dense or sparse format. - - Parameters - ---------- - out : Numpy.ndarray, optional - If given, the output will be written in-place into this array. - is_sparse : bool, optional - If set to True the method returns the matrix as a Scipy sparse matrix, if set to false - it returns the full matrix as a Numpy.ndarray - format : string, optional - Only relevant if is_sparse is True. Specifies the format in which the sparse matrix is to be stored. - Choose from "csr" (Compressed Sparse Row, default),"csc" (Compressed Sparse Column), "bsr" (Block Sparse Row ), - "lil" (List of Lists), "dok" (Dictionary of Keys), "coo" (COOrdinate format) and "dia" (DIAgonal). + Transforms the linear operator into a matrix, which is stored in sparse csr format. Returns ------- out : Numpy.ndarray or scipy.sparse.csr.csr_matrix - The matrix form of the linear operator. If ran in parallel each rank gets the full + The matrix form of the linear operator. If ran in parallel each rank gets the local matrix representation of the linear operator. """ # v will be the unit vector with which we compute Av = ith column of A. @@ -252,8 +240,7 @@ def tosparse(self): tmp2 = self.codomain.zeros() #We need to determine if we are a blockvector or a stencilvector but we are not able to use - #the BlockVectorSpace and StencilVectorSpace classes in here. So we use the following - #char. + #the BlockVectorSpace and StencilVectorSpace classes in here. So we use the following char. #This char is either a b for Blockvectors or an s for Stencilvectors BoS= str(self.domain)[15] @@ -412,26 +399,19 @@ def tosparse(self): return sparse.csr_matrix((data, (row, colarr)), shape=(numrows, numcols)) # Function that returns the matrix corresponding to the linear operator. Returns a numpy array. - def toarray(self, out=None, is_sparse=False, format="csr"): + def toarray(self, out=None): """ - Transforms the linear operator into a matrix, which is either stored in dense or sparse format. + Transforms the linear operator into a matrix, which is stored in dense format. Parameters ---------- out : Numpy.ndarray, optional If given, the output will be written in-place into this array. - is_sparse : bool, optional - If set to True the method returns the matrix as a Scipy sparse matrix, if set to false - it returns the full matrix as a Numpy.ndarray - format : string, optional - Only relevant if is_sparse is True. Specifies the format in which the sparse matrix is to be stored. - Choose from "csr" (Compressed Sparse Row, default),"csc" (Compressed Sparse Column), "bsr" (Block Sparse Row ), - "lil" (List of Lists), "dok" (Dictionary of Keys), "coo" (COOrdinate format) and "dia" (DIAgonal). - + Returns ------- - out : Numpy.ndarray or scipy.sparse.csr.csr_matrix - The matrix form of the linear operator. If ran in parallel each rank gets the full + out : Numpy.ndarray + The matrix form of the linear operator. If ran in parallel each rank gets the local matrix representation of the linear operator. """ # v will be the unit vector with which we compute Av = ith column of A. @@ -451,26 +431,15 @@ def toarray(self, out=None, is_sparse=False, format="csr"): rank = comm.Get_rank() size = comm.Get_size() - if (is_sparse == False): - if out is None: - # We declare the matrix form of our linear operator - out = np.zeros( - [self.codomain.dimension, self.domain.dimension], dtype=self.dtype) - else: - assert isinstance(out, np.ndarray) - assert out.shape[0] == self.codomain.dimension - assert out.shape[1] == self.domain.dimension + if out is None: + # We declare the matrix form of our linear operator + out = np.zeros( + [self.codomain.dimension, self.domain.dimension], dtype=self.dtype) else: - if out is not None: - raise Exception( - 'If is_sparse is True then out must be set to None.') - numrows = self.codomain.dimension - numcols = self.domain.dimension - # We define a list to store the non-zero data, a list to sotre the row index of said data and a list to store the column index. - data = [] - row = [] - colarr = [] - + assert isinstance(out, np.ndarray) + assert out.shape[0] == self.codomain.dimension + assert out.shape[1] == self.domain.dimension + # V is either a BlockVector or a StencilVector depending on the domain of the linear operator. if BoS == "b": # we collect all starts and ends in two big lists @@ -531,15 +500,7 @@ def toarray(self, out=None, is_sparse=False, format="csr"): # Compute to which column this iteration belongs col = spoint col += np.ravel_multi_index(i, npts[h]) - if is_sparse == False: - out[:, col] = tmp2.toarray() - else: - aux = tmp2.toarray() - # We now need to now which entries on tmp2 are non-zero and store then in our data list - for l in np.where(aux != 0)[0]: - data.append(aux[l]) - colarr.append(col) - row.append(l) + out[:, col] = tmp2.toarray() if (rank == currentrank): v[h][i] = 0.0 v[h].update_ghost_regions() @@ -597,15 +558,7 @@ def toarray(self, out=None, is_sparse=False, format="csr"): self.dot(v, out=tmp2) # Compute to which column this iteration belongs col = np.ravel_multi_index(i, npts) - if is_sparse == False: - out[:, col] = tmp2.toarray() - else: - aux = tmp2.toarray() - # We now need to now which entries on tmp2 are non-zero and store then in our data list - for l in np.where(aux != 0)[0]: - data.append(aux[l]) - colarr.append(col) - row.append(l) + out[:, col] = tmp2.toarray() if (rank == currentrank): v[i] = 0.0 v.update_ghost_regions() @@ -614,59 +567,8 @@ def toarray(self, out=None, is_sparse=False, format="csr"): # I cannot conceive any situation where this error should be thrown, but I put it here just in case something unexpected happens. raise Exception( 'Function toarray_struphy() only supports Stencil Vectors or Block Vectors.') - - if is_sparse == False: - return out - else: - # Gather all rows on rank 0 - gathered_rows = comm.gather(row, root=0) - # Gather all colarr on rank 0 - gathered_cols = comm.gather(colarr, root=0) - # Gather all data on rank 0 - gathered_data = comm.gather(data, root=0) - - if rank == 0: - # Rank 0 collects all rows from other ranks - all_rows = [ - item for sublist in gathered_rows for item in sublist] - # Rank 0 collects all columns from other ranks - all_cols = [ - item for sublist in gathered_cols for item in sublist] - # Rank 0 collects all data from other ranks - all_data = [ - item for sublist in gathered_data for item in sublist] - - # Broadcast 'all_rows' to all other ranks - comm.bcast(all_rows, root=0) - # Broadcast 'all_cols' to all other ranks - comm.bcast(all_cols, root=0) - # Broadcast 'all_data' to all other ranks - comm.bcast(all_data, root=0) - else: - # Other ranks receive the 'all_rows' list through broadcast - all_rows = comm.bcast(None, root=0) - # Other ranks receive the 'all_cols' list through broadcast - all_cols = comm.bcast(None, root=0) - # Other ranks receive the 'all_data' list through broadcast - all_data = comm.bcast(None, root=0) - - if format == "csr": - return sparse.csr_matrix((all_data, (all_rows, all_cols)), shape=(numrows, numcols)) - elif format == "csc": - return sparse.csc_matrix((all_data, (all_rows, all_cols)), shape=(numrows, numcols)) - elif format == "bsr": - return sparse.bsr_matrix((all_data, (all_rows, all_cols)), shape=(numrows, numcols)) - elif format == "lil": - return sparse.csr_matrix((all_data, (all_rows, all_cols)), shape=(numrows, numcols)).tolil() - elif format == "dok": - return sparse.csr_matrix((all_data, (all_rows, all_cols)), shape=(numrows, numcols)).todok() - elif format == "coo": - return sparse.coo_matrix((all_data, (all_rows, all_cols)), shape=(numrows, numcols)) - elif format == "dia": - return sparse.csr_matrix((all_data, (all_rows, all_cols)), shape=(numrows, numcols)).todia() - else: - raise Exception( - 'The selected sparse matrix format must be one of the following : csr, csc, bsr, lil, dok, coo or dia.') + return out + @abstractmethod diff --git a/psydac/linalg/tests/test_toarray.py b/psydac/linalg/tests/test_toarray.py index f347613ce..db29c0257 100644 --- a/psydac/linalg/tests/test_toarray.py +++ b/psydac/linalg/tests/test_toarray.py @@ -52,10 +52,6 @@ def get_StencilVectorSpace(n1, n2, p1, p2, P1, P2): def test_square_stencil_basic(n1, n2, p1, p2, P1=False, P2=False): - # 1. Initiate square LOs S,S1 (StencilMatrix), I (IdentityOperator), Z (ZeroOperator) and a Stencilvector v - # 2. Test general basic operations - # 3. Test special cases - ### ### 1. Initiation ### @@ -116,6 +112,7 @@ def test_square_stencil_basic(n1, n2, p1, p2, P1=False, P2=False): Apoarr = np.matmul(A1, A1) Apoarr = np.matmul(Apoarr, A1) + assert isinstance(Spo3, PowerLinearOperator) assert np.array_equal(Spo3arr, Apoarr) ### @@ -140,11 +137,6 @@ def test_square_stencil_basic(n1, n2, p1, p2, P1=False, P2=False): def test_square_block_basic(n1, n2, p1, p2, P1=False, P2=False): - # 1. Initiate square LOs S,S1 (StencilMatrix), Z (ZeroOperator) and a Stencilvector v - # Initiate square LOs B,B1 (BlockLO), BZ (ZeroOperator), BI (IdentityOperator) and a BlockVector vb - # 2. Test general basic operations - # 3. Test special cases - # Initiate StencilVectorSpace V = get_StencilVectorSpace(n1, n2, p1, p2, P1, P2) diff --git a/psydac/linalg/tests/test_toarray_parallel.py b/psydac/linalg/tests/test_toarray_parallel.py index 9580f3cfe..daaae2c25 100644 --- a/psydac/linalg/tests/test_toarray_parallel.py +++ b/psydac/linalg/tests/test_toarray_parallel.py @@ -155,8 +155,7 @@ def test_block_linear_operator_parallel_dot( dtype, n1, n2, p1, p2, P1, P2 ): aux = i.flatten() for j in aux: v1arr.append(j) - - # Test copy with an out + # Create random matrix N1 = StencilMatrix( V, V ) N2 = StencilMatrix( V, V ) @@ -179,7 +178,9 @@ def test_block_linear_operator_parallel_dot( dtype, n1, n2, p1, p2, P1, P2 ): #Now we use a PowerLinearOperator KP = K**3 + #We compute the local matrix representation of KP KParrloc = KP.toarray() + #We compute the global matrix representation of KP KParr = np.zeros(np.shape(KParrloc)) comm.Allreduce(KParrloc, KParr, op=MPI.SUM) diff --git a/psydac/linalg/tests/test_tosparse_parallel.py b/psydac/linalg/tests/test_tosparse_parallel.py index 4eb4a5c46..6c574bed8 100644 --- a/psydac/linalg/tests/test_tosparse_parallel.py +++ b/psydac/linalg/tests/test_tosparse_parallel.py @@ -156,8 +156,7 @@ def test_block_linear_operator_parallel_dot( dtype, n1, n2, p1, p2, P1, P2 ): aux = i.flatten() for j in aux: v1arr.append(j) - - # Test copy with an out + # Create random matrix N1 = StencilMatrix( V, V ) N2 = StencilMatrix( V, V ) From 564005f3e10c26b42aaa7fb2f53a3a2b5d84f55c Mon Sep 17 00:00:00 2001 From: mateomarin97 Date: Tue, 12 Dec 2023 16:55:22 +0000 Subject: [PATCH 03/18] Implemented some of the suggested changes --- psydac/linalg/basic.py | 69 +++++++++++-------- psydac/linalg/tests/test_toarray_parallel.py | 53 ++++++++++++-- psydac/linalg/tests/test_tosparse_parallel.py | 56 ++++++++++++--- 3 files changed, 134 insertions(+), 44 deletions(-) diff --git a/psydac/linalg/basic.py b/psydac/linalg/basic.py index e31d975a2..ff521375c 100644 --- a/psydac/linalg/basic.py +++ b/psydac/linalg/basic.py @@ -240,14 +240,23 @@ def tosparse(self): tmp2 = self.codomain.zeros() #We need to determine if we are a blockvector or a stencilvector but we are not able to use - #the BlockVectorSpace and StencilVectorSpace classes in here. So we use the following char. - #This char is either a b for Blockvectors or an s for Stencilvectors - BoS= str(self.domain)[15] - + #the BlockVectorSpace and StencilVectorSpace classes in here. So we check if domain has the spaces + #attribute in which case the domain would be a BlockVectorSpace. If that is not the case we check + #if the domain has the cart atrribute, in which case it will be a StencilVectorSpace. + + if hasattr(self.domain, 'spaces'): + BoS = "b" + elif hasattr(self.domain, 'cart'): + BoS = "s" + else: + raise Exception( + 'The domain of the LinearOperator must be a BlockVectorSpace or a StencilVectorSpace.') + if BoS == "b": comm = self.domain.spaces[0].cart.comm elif BoS == "s": comm = self.domain.cart.comm + rank = comm.Get_rank() size = comm.Get_size() @@ -303,12 +312,12 @@ def tosparse(self): npredim = 0 # We iterate over the stencil vectors inside the BlockVector for h in range(nsp): - itterables = [] + iterables = [] for i in range(ndim[h]): - itterables.append( + iterables.append( range(allstarts[currentrank][i+npredim], allends[currentrank][i+npredim]+1)) # We iterate over all the entries that belong to rank number currentrank - for i in itertools.product(*itterables): + for i in itertools.product(*iterables): if (rank == currentrank): v[h][i] = 1.0 v[h].update_ghost_regions() @@ -368,12 +377,12 @@ def tosparse(self): currentrank = 0 # Each rank will take care of setting to 1 each one of its entries while all other entries remain zero. while (currentrank < size): - itterables = [] + iterables = [] for i in range(ndim): - itterables.append( + iterables.append( range(allstarts[currentrank][i], allends[currentrank][i]+1)) # We iterate over all the entries that belong to rank number currentrank - for i in itertools.product(*itterables): + for i in itertools.product(*iterables): if (rank == currentrank): v[i] = 1.0 v.update_ghost_regions() @@ -391,10 +400,6 @@ def tosparse(self): v[i] = 0.0 v.update_ghost_regions() currentrank += 1 - else: - # I cannot conceive any situation where this error should be thrown, but I put it here just in case something unexpected happens. - raise Exception( - 'Function toarray_struphy() only supports Stencil Vectors or Block Vectors.') return sparse.csr_matrix((data, (row, colarr)), shape=(numrows, numcols)) @@ -418,16 +423,25 @@ def toarray(self, out=None): v = self.domain.zeros() # We define a temporal vector tmp2 = self.codomain.zeros() + #We need to determine if we are a blockvector or a stencilvector but we are not able to use - #the BlockVectorSpace and StencilVectorSpace classes in here. So we use the following - #char. - #This char is either a b for Blockvectors or an s for Stencilvectors - BoS= str(self.domain)[15] - + #the BlockVectorSpace and StencilVectorSpace classes in here. So we check if domain has the spaces + #attribute in which case the domain would be a BlockVectorSpace. If that is not the case we check + #if the domain has the cart atrribute, in which case it will be a StencilVectorSpace. + + if hasattr(self.domain, 'spaces'): + BoS = "b" + elif hasattr(self.domain, 'cart'): + BoS = "s" + else: + raise Exception( + 'The domain of the LinearOperator must be a BlockVectorSpace or a StencilVectorSpace.') + if BoS == "b": comm = self.domain.spaces[0].cart.comm elif BoS == "s": comm = self.domain.cart.comm + rank = comm.Get_rank() size = comm.Get_size() @@ -485,12 +499,12 @@ def toarray(self, out=None): npredim = 0 # We iterate over the stencil vectors inside the BlockVector for h in range(nsp): - itterables = [] + iterables = [] for i in range(ndim[h]): - itterables.append( + iterables.append( range(allstarts[currentrank][i+npredim], allends[currentrank][i+npredim]+1)) # We iterate over all the entries that belong to rank number currentrank - for i in itertools.product(*itterables): + for i in itertools.product(*iterables): if (rank == currentrank): v[h][i] = 1.0 v[h].update_ghost_regions() @@ -545,12 +559,12 @@ def toarray(self, out=None): currentrank = 0 # Each rank will take care of setting to 1 each one of its entries while all other entries remain zero. while (currentrank < size): - itterables = [] + iterables = [] for i in range(ndim): - itterables.append( + iterables.append( range(allstarts[currentrank][i], allends[currentrank][i]+1)) # We iterate over all the entries that belong to rank number currentrank - for i in itertools.product(*itterables): + for i in itertools.product(*iterables): if (rank == currentrank): v[i] = 1.0 v.update_ghost_regions() @@ -563,10 +577,7 @@ def toarray(self, out=None): v[i] = 0.0 v.update_ghost_regions() currentrank += 1 - else: - # I cannot conceive any situation where this error should be thrown, but I put it here just in case something unexpected happens. - raise Exception( - 'Function toarray_struphy() only supports Stencil Vectors or Block Vectors.') + return out diff --git a/psydac/linalg/tests/test_toarray_parallel.py b/psydac/linalg/tests/test_toarray_parallel.py index daaae2c25..8e24f8cf2 100644 --- a/psydac/linalg/tests/test_toarray_parallel.py +++ b/psydac/linalg/tests/test_toarray_parallel.py @@ -93,19 +93,35 @@ def compute_global_starts_ends(domain_decomposition, npts): def create_equal_random_arrays(W, seedv =123): np.random.seed(seedv) arr = [] - arr_psy = BlockVector(W) + if isinstance(W, BlockVectorSpace): + arr_psy = BlockVector(W) - for d, block in enumerate(arr_psy.blocks): + for d, block in enumerate(arr_psy.blocks): - dims = W.spaces[d].npts + dims = W.spaces[d].npts + + arr += [np.random.rand(*dims)] + + s = block.starts + e = block.ends + + arr_psy[d][s[0]:e[0] + 1, s[1]:e[1] + 1] = arr[-1][s[0]:e[0] + 1, s[1]:e[1] + 1] + arr_psy.update_ghost_regions() + elif isinstance(W, StencilVectorSpace): + arr_psy = StencilVector(W) + + dims = W.npts arr += [np.random.rand(*dims)] - s = block.starts - e = block.ends + s = arr_psy.starts + e = arr_psy.ends - arr_psy[d][s[0]:e[0] + 1, s[1]:e[1] + 1] = arr[-1][s[0]:e[0] + 1, s[1]:e[1] + 1] - arr_psy.update_ghost_regions() + arr_psy[s[0]:e[0] + 1, s[1]:e[1] + 1] = arr[-1][s[0]:e[0] + 1, s[1]:e[1] + 1] + + else: + raise Exception( + 'W must be a BlockVectorSpace or a StencilVectorSpace.') return arr, arr_psy @@ -146,6 +162,11 @@ def test_block_linear_operator_parallel_dot( dtype, n1, n2, p1, p2, P1, P2 ): # Create vector space, stencil matrix, and stencil vector V = StencilVectorSpace( cart, dtype=dtype ) + + v0arr , v0 = create_equal_random_arrays(V, seedv=4568) + + v0arr = v0arr[0].flatten() + # Create and Fill Block objects W = BlockVectorSpace(V, V) @@ -176,6 +197,24 @@ def test_block_linear_operator_parallel_dot( dtype, n1, n2, p1, p2, P1, P2 ): K[1,0] = N3 K[1,1] = N4 + ##### + #Testing toarray() for LinearOperators with domain being StencilVectorSpace + + #We take a power of the Stencil Matrix N1 + N1P = N1**3 + #We compute the local matrix representation of N1P + N1arrloc = N1P.toarray() + #We get the global matrix of the PowerLinearOperator + N1arr = np.zeros(np.shape(N1arrloc)) + comm.Allreduce(N1arrloc, N1arr, op=MPI.SUM) + + assert isinstance(N1P, PowerLinearOperator) + assert isinstance(N1P.domain, StencilVectorSpace) + compare_arrays(N1P.dot(v0), np.matmul(N1arr,v0arr), rank) + + ##### + #Testing toarray() for LinearOperators with domain being BlockVectorSpace + #Now we use a PowerLinearOperator KP = K**3 #We compute the local matrix representation of KP diff --git a/psydac/linalg/tests/test_tosparse_parallel.py b/psydac/linalg/tests/test_tosparse_parallel.py index 6c574bed8..f726c6259 100644 --- a/psydac/linalg/tests/test_tosparse_parallel.py +++ b/psydac/linalg/tests/test_tosparse_parallel.py @@ -93,19 +93,35 @@ def compute_global_starts_ends(domain_decomposition, npts): def create_equal_random_arrays(W, seedv =123): np.random.seed(seedv) arr = [] - arr_psy = BlockVector(W) + if isinstance(W, BlockVectorSpace): + arr_psy = BlockVector(W) - for d, block in enumerate(arr_psy.blocks): + for d, block in enumerate(arr_psy.blocks): - dims = W.spaces[d].npts + dims = W.spaces[d].npts + + arr += [np.random.rand(*dims)] + + s = block.starts + e = block.ends + + arr_psy[d][s[0]:e[0] + 1, s[1]:e[1] + 1] = arr[-1][s[0]:e[0] + 1, s[1]:e[1] + 1] + arr_psy.update_ghost_regions() + elif isinstance(W, StencilVectorSpace): + arr_psy = StencilVector(W) + + dims = W.npts arr += [np.random.rand(*dims)] - s = block.starts - e = block.ends + s = arr_psy.starts + e = arr_psy.ends - arr_psy[d][s[0]:e[0] + 1, s[1]:e[1] + 1] = arr[-1][s[0]:e[0] + 1, s[1]:e[1] + 1] - arr_psy.update_ghost_regions() + arr_psy[s[0]:e[0] + 1, s[1]:e[1] + 1] = arr[-1][s[0]:e[0] + 1, s[1]:e[1] + 1] + + else: + raise Exception( + 'W must be a BlockVectorSpace or a StencilVectorSpace.') return arr, arr_psy @@ -147,6 +163,10 @@ def test_block_linear_operator_parallel_dot( dtype, n1, n2, p1, p2, P1, P2 ): # Create vector space, stencil matrix, and stencil vector V = StencilVectorSpace( cart, dtype=dtype ) + v0arr , v0 = create_equal_random_arrays(V, seedv=4568) + + v0arr = v0arr[0].flatten() + # Create and Fill Block objects W = BlockVectorSpace(V, V) @@ -177,7 +197,26 @@ def test_block_linear_operator_parallel_dot( dtype, n1, n2, p1, p2, P1, P2 ): K[1,0] = N3 K[1,1] = N4 - #We take a power of the LinearOperator + + ##### + #Testing tosparse() for LinearOperators with domain being StencilVectorSpace + + #We take a power of the Stencil Matrix N1 + N1P = N1**3 + #We turn our PowerLinearOperator into a sparse matrix, and then into a numpy array + N1arrloc = N1P.tosparse().toarray() + #We get the global matrix of the PowerLinearOperator + N1arr = np.zeros(np.shape(N1arrloc)) + comm.Allreduce(N1arrloc, N1arr, op=MPI.SUM) + + assert isinstance(N1P, PowerLinearOperator) + assert isinstance(N1P.domain, StencilVectorSpace) + compare_arrays(N1P.dot(v0), N1arr.dot(v0arr), rank) + + ##### + #Testing tosparse() for LinearOperators with domain being BlockVectorSpace + + #We take a power of the BlockLinearOperator KP = K**3 #We turn our PowerLinearOperator into a sparse matrix, and then into a numpy array KLarrloc = KP.tosparse().toarray() @@ -186,6 +225,7 @@ def test_block_linear_operator_parallel_dot( dtype, n1, n2, p1, p2, P1, P2 ): comm.Allreduce(KLarrloc, KLarr, op=MPI.SUM) assert isinstance(KP, PowerLinearOperator) + assert isinstance(KP.domain, BlockVectorSpace) compare_arrays(KP.dot(v1), KLarr.dot(v1arr), rank) From f5de4d3b625861503419f877908516460f6f5a9e Mon Sep 17 00:00:00 2001 From: mateomarin97 Date: Wed, 13 Dec 2023 15:50:10 +0000 Subject: [PATCH 04/18] Reduced the amount of code by writing and additional private function _tosparse_array, that contains the functionality to get both the sparse matrix and the full matrix --- psydac/linalg/basic.py | 421 ++++++------------ psydac/linalg/tests/test_toarray_parallel.py | 10 +- psydac/linalg/tests/test_tosparse_parallel.py | 12 +- 3 files changed, 147 insertions(+), 296 deletions(-) diff --git a/psydac/linalg/basic.py b/psydac/linalg/basic.py index ff521375c..7d7f90053 100644 --- a/psydac/linalg/basic.py +++ b/psydac/linalg/basic.py @@ -222,28 +222,35 @@ def codomain(self): @abstractmethod def dtype(self): pass - - # Function that returns the local matrix corresponding to the linear operator. Returns a scipy.sparse.csr.csr_matrix. - def tosparse(self): + + #Private function that contains the functionality to transform a linear operator into a scipy.sparse.csr.csr_matrix or a numpy array + def __tosparse_array(self, out=None, is_sparse=False): """ - Transforms the linear operator into a matrix, which is stored in sparse csr format. + Transforms the linear operator into a matrix, which is either stored in dense or sparse format. + + Parameters + ---------- + out : Numpy.ndarray, optional + If given, the output will be written in-place into this array. + is_sparse : bool, optional + If set to True the method returns the matrix as a Scipy sparse matrix, if set to false + it returns the full matrix as a Numpy.ndarray Returns ------- out : Numpy.ndarray or scipy.sparse.csr.csr_matrix - The matrix form of the linear operator. If ran in parallel each rank gets the local + The matrix form of the linear operator. If ran in parallel each rank gets the full matrix representation of the linear operator. """ # v will be the unit vector with which we compute Av = ith column of A. v = self.domain.zeros() # We define a temporal vector tmp2 = self.codomain.zeros() - + #We need to determine if we are a blockvector or a stencilvector but we are not able to use #the BlockVectorSpace and StencilVectorSpace classes in here. So we check if domain has the spaces #attribute in which case the domain would be a BlockVectorSpace. If that is not the case we check #if the domain has the cart atrribute, in which case it will be a StencilVectorSpace. - if hasattr(self.domain, 'spaces'): BoS = "b" elif hasattr(self.domain, 'cart'): @@ -256,16 +263,28 @@ def tosparse(self): comm = self.domain.spaces[0].cart.comm elif BoS == "s": comm = self.domain.cart.comm - rank = comm.Get_rank() size = comm.Get_size() - - numrows = self.codomain.dimension - numcols = self.domain.dimension - # We define a list to store the non-zero data, a list to sotre the row index of said data and a list to store the column index. - data = [] - row = [] - colarr = [] + + if (is_sparse == False): + if out is None: + # We declare the matrix form of our linear operator + out = np.zeros( + [self.codomain.dimension, self.domain.dimension], dtype=self.dtype) + else: + assert isinstance(out, np.ndarray) + assert out.shape[0] == self.codomain.dimension + assert out.shape[1] == self.domain.dimension + else: + if out is not None: + raise Exception( + 'If is_sparse is True then out must be set to None.') + numrows = self.codomain.dimension + numcols = self.domain.dimension + # We define a list to store the non-zero data, a list to sotre the row index of said data and a list to store the column index. + data = [] + row = [] + colarr = [] # V is either a BlockVector or a StencilVector depending on the domain of the linear operator. if BoS == "b": @@ -278,131 +297,121 @@ def tosparse(self): nsp = len(self.domain.spaces) # We get the number of dimensions each space has. ndim = [sp.ndim for sp in self.domain.spaces] - - # First each rank is going to need to know the starts and ends of all other ranks - startsarr = np.array([starts[i][j] for i in range(nsp) - for j in range(ndim[i])], dtype=int) - - # Create an array to store gathered data from all ranks - allstarts = np.empty(size * len(startsarr), dtype=int) - - # Use Allgather to gather 'starts' from all ranks into 'allstarts' - comm.Allgather(startsarr, allstarts) - - # Reshape 'allstarts' to have 9 columns and 'size' rows - allstarts = allstarts.reshape((size, len(startsarr))) - - endsarr = np.array([ends[i][j] for i in range(nsp) - for j in range(ndim[i])], dtype=int) - # Create an array to store gathered data from all ranks - allends = np.empty(size * len(endsarr), dtype=int) - - # Use Allgather to gather 'ends' from all ranks into 'allends' - comm.Allgather(endsarr, allends) - - # Reshape 'allends' to have 9 columns and 'size' rows - allends = allends.reshape((size, len(endsarr))) - - currentrank = 0 - # Each rank will take care of setting to 1 each one of its entries while all other entries remain zero. - while (currentrank < size): - # since the size of npts changes denpending on h we need to compute a starting point for - # our column index - spoint = 0 - npredim = 0 - # We iterate over the stencil vectors inside the BlockVector - for h in range(nsp): - iterables = [] - for i in range(ndim[h]): - iterables.append( - range(allstarts[currentrank][i+npredim], allends[currentrank][i+npredim]+1)) - # We iterate over all the entries that belong to rank number currentrank - for i in itertools.product(*iterables): + elif BoS == "s": + # We get the start and endpoint for each sublist in v + starts = [v.starts] + ends = [v.ends] + # We get the dimensions of the StencilVector + npts = [self.domain.npts] + # We get the number of space we have + nsp = 1 + # We get the number of dimensions the StencilVectorSpace has. + ndim = [self.domain.ndim] + + # First each rank is going to need to know the starts and ends of all other ranks + startsarr = np.array([starts[i][j] for i in range(nsp) + for j in range(ndim[i])], dtype=int) + + endsarr = np.array([ends[i][j] for i in range(nsp) + for j in range(ndim[i])], dtype=int) + + # Create an array to store gathered data from all ranks + allstarts = np.empty(size * len(startsarr), dtype=int) + + # Use Allgather to gather 'starts' from all ranks into 'allstarts' + comm.Allgather(startsarr, allstarts) + + # Reshape 'allstarts' to have 9 columns and 'size' rows + allstarts = allstarts.reshape((size, len(startsarr))) + + # Create an array to store gathered data from all ranks + allends = np.empty(size * len(endsarr), dtype=int) + + # Use Allgather to gather 'ends' from all ranks into 'allends' + comm.Allgather(endsarr, allends) + + # Reshape 'allends' to have 9 columns and 'size' rows + allends = allends.reshape((size, len(endsarr))) + + currentrank = 0 + # Each rank will take care of setting to 1 each one of its entries while all other entries remain zero. + while (currentrank < size): + # since the size of npts changes denpending on h we need to compute a starting point for + # our column index + spoint = 0 + npredim = 0 + # We iterate over the stencil vectors inside the BlockVector + for h in range(nsp): + itterables = [] + for i in range(ndim[h]): + itterables.append( + range(allstarts[currentrank][i+npredim], allends[currentrank][i+npredim]+1)) + # We iterate over all the entries that belong to rank number currentrank + for i in itertools.product(*itterables): + + ######################################### + if BoS == "b": if (rank == currentrank): v[h][i] = 1.0 v[h].update_ghost_regions() - # Compute dot product with the linear operator. - tmp2 *= 0. - self.dot(v, out=tmp2) - # Compute to which column this iteration belongs - col = spoint - col += np.ravel_multi_index(i, npts[h]) + elif BoS == "s": + if (rank == currentrank): + v[i] = 1.0 + v.update_ghost_regions() + ######################################### + + # Compute dot product with the linear operator. + self.dot(v, out=tmp2) + # Compute to which column this iteration belongs + col = spoint + col += np.ravel_multi_index(i, npts[h]) + if is_sparse == False: + out[:, col] = tmp2.toarray() + else: aux = tmp2.toarray() # We now need to now which entries on tmp2 are non-zero and store then in our data list for l in np.where(aux != 0)[0]: data.append(aux[l]) colarr.append(col) row.append(l) + + ################################# + if BoS == "b": if (rank == currentrank): v[h][i] = 0.0 v[h].update_ghost_regions() - cummulative = 1 - for i in range(ndim[h]): - cummulative *= npts[h][i] - spoint += cummulative - npredim += ndim[h] - currentrank += 1 - elif BoS == "s": - # We get the start and endpoint for each sublist in v - starts = v.starts - ends = v.ends - # We get the dimensions of the StencilVector - npts = self.domain.npts - # We get the number of space we have - nsp = 1 - # We get the number of dimensions the StencilVectorSpace has. - ndim = self.domain.ndim - - # First each rank is going to need to know the starts and ends of all other ranks - startsarr = np.array([starts[j] for j in range(ndim)], dtype=int) - # Create an array to store gathered data from all ranks - allstarts = np.empty(size * len(startsarr), dtype=int) - - # Use Allgather to gather 'starts' from all ranks into 'allstarts' - comm.Allgather(startsarr, allstarts) - - # Reshape 'allstarts' to have 3 columns and 'size' rows - allstarts = allstarts.reshape((size, len(startsarr))) - - endsarr = np.array([ends[j] for j in range(ndim)], dtype=int) - # Create an array to store gathered data from all ranks - allends = np.empty(size * len(endsarr), dtype=int) - - # Use Allgather to gather 'ends' from all ranks into 'allends' - comm.Allgather(endsarr, allends) - - # Reshape 'allends' to have 3 columns and 'size' rows - allends = allends.reshape((size, len(endsarr))) - - currentrank = 0 - # Each rank will take care of setting to 1 each one of its entries while all other entries remain zero. - while (currentrank < size): - iterables = [] - for i in range(ndim): - iterables.append( - range(allstarts[currentrank][i], allends[currentrank][i]+1)) - # We iterate over all the entries that belong to rank number currentrank - for i in itertools.product(*iterables): - if (rank == currentrank): - v[i] = 1.0 - v.update_ghost_regions() - # Compute dot product with the linear operator. - self.dot(v, out=tmp2) - # Compute to which column this iteration belongs - col = np.ravel_multi_index(i, npts) - aux = tmp2.toarray() - # We now need to now which entries on tmp2 are non-zero and store then in our data list - for l in np.where(aux != 0)[0]: - data.append(aux[l]) - colarr.append(col) - row.append(l) - if (rank == currentrank): - v[i] = 0.0 - v.update_ghost_regions() - currentrank += 1 + elif BoS == "s": + if (rank == currentrank): + v[i] = 0.0 + v.update_ghost_regions() + ################################## + cummulative = 1 + for i in range(ndim[h]): + cummulative *= npts[h][i] + spoint += cummulative + npredim += ndim[h] + currentrank += 1 - return sparse.csr_matrix((data, (row, colarr)), shape=(numrows, numcols)) + if is_sparse == False: + return out + else: + return sparse.csr_matrix((data, (row, colarr)), shape=(numrows, numcols)) + + + # Function that returns the local matrix corresponding to the linear operator. Returns a scipy.sparse.csr.csr_matrix. + def tosparse(self): + """ + Transforms the linear operator into a matrix, which is stored in sparse csr format. + Returns + ------- + out : Numpy.ndarray or scipy.sparse.csr.csr_matrix + The matrix form of the linear operator. If ran in parallel each rank gets the local + matrix representation of the linear operator. + """ + return self.__tosparse_array(is_sparse=True) + + # Function that returns the matrix corresponding to the linear operator. Returns a numpy array. def toarray(self, out=None): """ @@ -419,166 +428,8 @@ def toarray(self, out=None): The matrix form of the linear operator. If ran in parallel each rank gets the local matrix representation of the linear operator. """ - # v will be the unit vector with which we compute Av = ith column of A. - v = self.domain.zeros() - # We define a temporal vector - tmp2 = self.codomain.zeros() + return self.__tosparse_array(out=out, is_sparse=False) - #We need to determine if we are a blockvector or a stencilvector but we are not able to use - #the BlockVectorSpace and StencilVectorSpace classes in here. So we check if domain has the spaces - #attribute in which case the domain would be a BlockVectorSpace. If that is not the case we check - #if the domain has the cart atrribute, in which case it will be a StencilVectorSpace. - - if hasattr(self.domain, 'spaces'): - BoS = "b" - elif hasattr(self.domain, 'cart'): - BoS = "s" - else: - raise Exception( - 'The domain of the LinearOperator must be a BlockVectorSpace or a StencilVectorSpace.') - - if BoS == "b": - comm = self.domain.spaces[0].cart.comm - elif BoS == "s": - comm = self.domain.cart.comm - - rank = comm.Get_rank() - size = comm.Get_size() - - if out is None: - # We declare the matrix form of our linear operator - out = np.zeros( - [self.codomain.dimension, self.domain.dimension], dtype=self.dtype) - else: - assert isinstance(out, np.ndarray) - assert out.shape[0] == self.codomain.dimension - assert out.shape[1] == self.domain.dimension - - # V is either a BlockVector or a StencilVector depending on the domain of the linear operator. - if BoS == "b": - # we collect all starts and ends in two big lists - starts = [vi.starts for vi in v] - ends = [vi.ends for vi in v] - # We collect the dimension of the BlockVector - npts = [sp.npts for sp in self.domain.spaces] - # We get the number of space we have - nsp = len(self.domain.spaces) - # We get the number of dimensions each space has. - ndim = [sp.ndim for sp in self.domain.spaces] - - # First each rank is going to need to know the starts and ends of all other ranks - startsarr = np.array([starts[i][j] for i in range(nsp) - for j in range(ndim[i])], dtype=int) - - # Create an array to store gathered data from all ranks - allstarts = np.empty(size * len(startsarr), dtype=int) - - # Use Allgather to gather 'starts' from all ranks into 'allstarts' - comm.Allgather(startsarr, allstarts) - - # Reshape 'allstarts' to have 9 columns and 'size' rows - allstarts = allstarts.reshape((size, len(startsarr))) - - endsarr = np.array([ends[i][j] for i in range(nsp) - for j in range(ndim[i])], dtype=int) - # Create an array to store gathered data from all ranks - allends = np.empty(size * len(endsarr), dtype=int) - - # Use Allgather to gather 'ends' from all ranks into 'allends' - comm.Allgather(endsarr, allends) - - # Reshape 'allends' to have 9 columns and 'size' rows - allends = allends.reshape((size, len(endsarr))) - - currentrank = 0 - # Each rank will take care of setting to 1 each one of its entries while all other entries remain zero. - while (currentrank < size): - # since the size of npts changes denpending on h we need to compute a starting point for - # our column index - spoint = 0 - npredim = 0 - # We iterate over the stencil vectors inside the BlockVector - for h in range(nsp): - iterables = [] - for i in range(ndim[h]): - iterables.append( - range(allstarts[currentrank][i+npredim], allends[currentrank][i+npredim]+1)) - # We iterate over all the entries that belong to rank number currentrank - for i in itertools.product(*iterables): - if (rank == currentrank): - v[h][i] = 1.0 - v[h].update_ghost_regions() - # Compute dot product with the linear operator. - tmp2 *= 0. - self.dot(v, out=tmp2) - # Compute to which column this iteration belongs - col = spoint - col += np.ravel_multi_index(i, npts[h]) - out[:, col] = tmp2.toarray() - if (rank == currentrank): - v[h][i] = 0.0 - v[h].update_ghost_regions() - cummulative = 1 - for i in range(ndim[h]): - cummulative *= npts[h][i] - spoint += cummulative - npredim += ndim[h] - currentrank += 1 - elif BoS == "s": - # We get the start and endpoint for each sublist in v - starts = v.starts - ends = v.ends - # We get the dimensions of the StencilVector - npts = self.domain.npts - # We get the number of space we have - nsp = 1 - # We get the number of dimensions the StencilVectorSpace has. - ndim = self.domain.ndim - - # First each rank is going to need to know the starts and ends of all other ranks - startsarr = np.array([starts[j] for j in range(ndim)], dtype=int) - # Create an array to store gathered data from all ranks - allstarts = np.empty(size * len(startsarr), dtype=int) - - # Use Allgather to gather 'starts' from all ranks into 'allstarts' - comm.Allgather(startsarr, allstarts) - - # Reshape 'allstarts' to have 3 columns and 'size' rows - allstarts = allstarts.reshape((size, len(startsarr))) - - endsarr = np.array([ends[j] for j in range(ndim)], dtype=int) - # Create an array to store gathered data from all ranks - allends = np.empty(size * len(endsarr), dtype=int) - - # Use Allgather to gather 'ends' from all ranks into 'allends' - comm.Allgather(endsarr, allends) - - # Reshape 'allends' to have 3 columns and 'size' rows - allends = allends.reshape((size, len(endsarr))) - - currentrank = 0 - # Each rank will take care of setting to 1 each one of its entries while all other entries remain zero. - while (currentrank < size): - iterables = [] - for i in range(ndim): - iterables.append( - range(allstarts[currentrank][i], allends[currentrank][i]+1)) - # We iterate over all the entries that belong to rank number currentrank - for i in itertools.product(*iterables): - if (rank == currentrank): - v[i] = 1.0 - v.update_ghost_regions() - # Compute dot product with the linear operator. - self.dot(v, out=tmp2) - # Compute to which column this iteration belongs - col = np.ravel_multi_index(i, npts) - out[:, col] = tmp2.toarray() - if (rank == currentrank): - v[i] = 0.0 - v.update_ghost_regions() - currentrank += 1 - - return out diff --git a/psydac/linalg/tests/test_toarray_parallel.py b/psydac/linalg/tests/test_toarray_parallel.py index 8e24f8cf2..cbe4f77e9 100644 --- a/psydac/linalg/tests/test_toarray_parallel.py +++ b/psydac/linalg/tests/test_toarray_parallel.py @@ -128,9 +128,9 @@ def create_equal_random_arrays(W, seedv =123): @pytest.mark.parametrize( 'dtype', [float] ) -@pytest.mark.parametrize( 'n1', [8, 16] ) -@pytest.mark.parametrize( 'n2', [8, 32] ) -@pytest.mark.parametrize( 'p1', [1, 3] ) +@pytest.mark.parametrize( 'n1', [16, 32] ) +@pytest.mark.parametrize( 'n2', [16, 32] ) +@pytest.mark.parametrize( 'p1', [1, 2] ) @pytest.mark.parametrize( 'p2', [2] ) @pytest.mark.parametrize( 'P1', [True, False] ) @pytest.mark.parametrize( 'P2', [True] ) @@ -205,7 +205,7 @@ def test_block_linear_operator_parallel_dot( dtype, n1, n2, p1, p2, P1, P2 ): #We compute the local matrix representation of N1P N1arrloc = N1P.toarray() #We get the global matrix of the PowerLinearOperator - N1arr = np.zeros(np.shape(N1arrloc)) + N1arr = np.zeros(np.shape(N1arrloc), dtype = dtype) comm.Allreduce(N1arrloc, N1arr, op=MPI.SUM) assert isinstance(N1P, PowerLinearOperator) @@ -220,7 +220,7 @@ def test_block_linear_operator_parallel_dot( dtype, n1, n2, p1, p2, P1, P2 ): #We compute the local matrix representation of KP KParrloc = KP.toarray() #We compute the global matrix representation of KP - KParr = np.zeros(np.shape(KParrloc)) + KParr = np.zeros(np.shape(KParrloc), dtype = dtype) comm.Allreduce(KParrloc, KParr, op=MPI.SUM) diff --git a/psydac/linalg/tests/test_tosparse_parallel.py b/psydac/linalg/tests/test_tosparse_parallel.py index f726c6259..4ea9caf49 100644 --- a/psydac/linalg/tests/test_tosparse_parallel.py +++ b/psydac/linalg/tests/test_tosparse_parallel.py @@ -127,10 +127,10 @@ def create_equal_random_arrays(W, seedv =123): -@pytest.mark.parametrize( 'dtype', [float] ) -@pytest.mark.parametrize( 'n1', [8, 16] ) -@pytest.mark.parametrize( 'n2', [8, 32] ) -@pytest.mark.parametrize( 'p1', [1, 3] ) +@pytest.mark.parametrize( 'dtype', [float, complex] ) +@pytest.mark.parametrize( 'n1', [16] ) +@pytest.mark.parametrize( 'n2', [16, 32] ) +@pytest.mark.parametrize( 'p1', [1, 2] ) @pytest.mark.parametrize( 'p2', [2] ) @pytest.mark.parametrize( 'P1', [True, False] ) @pytest.mark.parametrize( 'P2', [True] ) @@ -206,7 +206,7 @@ def test_block_linear_operator_parallel_dot( dtype, n1, n2, p1, p2, P1, P2 ): #We turn our PowerLinearOperator into a sparse matrix, and then into a numpy array N1arrloc = N1P.tosparse().toarray() #We get the global matrix of the PowerLinearOperator - N1arr = np.zeros(np.shape(N1arrloc)) + N1arr = np.zeros(np.shape(N1arrloc), dtype=dtype) comm.Allreduce(N1arrloc, N1arr, op=MPI.SUM) assert isinstance(N1P, PowerLinearOperator) @@ -221,7 +221,7 @@ def test_block_linear_operator_parallel_dot( dtype, n1, n2, p1, p2, P1, P2 ): #We turn our PowerLinearOperator into a sparse matrix, and then into a numpy array KLarrloc = KP.tosparse().toarray() #We get the global matrix of the operator - KLarr = np.zeros(np.shape(KLarrloc)) + KLarr = np.zeros(np.shape(KLarrloc), dtype=dtype) comm.Allreduce(KLarrloc, KLarr, op=MPI.SUM) assert isinstance(KP, PowerLinearOperator) From 484365ad19c7a1996abdfa4bf3793bee9a83b5ee Mon Sep 17 00:00:00 2001 From: mateomarin97 Date: Thu, 18 Jan 2024 14:43:57 +0000 Subject: [PATCH 05/18] Optimized the toarray_sparse function by removing toarray() calls inside the main for loop that kept allocating memory unecessarily --- psydac/linalg/basic.py | 103 ++++++++++++++++++++++++++++++++++++----- 1 file changed, 91 insertions(+), 12 deletions(-) diff --git a/psydac/linalg/basic.py b/psydac/linalg/basic.py index 7d7f90053..163b4c3d4 100644 --- a/psydac/linalg/basic.py +++ b/psydac/linalg/basic.py @@ -223,7 +223,6 @@ def codomain(self): def dtype(self): pass - #Private function that contains the functionality to transform a linear operator into a scipy.sparse.csr.csr_matrix or a numpy array def __tosparse_array(self, out=None, is_sparse=False): """ Transforms the linear operator into a matrix, which is either stored in dense or sparse format. @@ -258,6 +257,15 @@ def __tosparse_array(self, out=None, is_sparse=False): else: raise Exception( 'The domain of the LinearOperator must be a BlockVectorSpace or a StencilVectorSpace.') + + #We also need to know if the codomain is a StencilVectorSpace or a BlockVectorSpace + if hasattr(self.codomain, 'spaces'): + BoS2 = "b" + elif hasattr(self.codomain, 'cart'): + BoS2 = "s" + else: + raise Exception( + 'The codomain of the LinearOperator must be a BlockVectorSpace or a StencilVectorSpace.') if BoS == "b": comm = self.domain.spaces[0].cart.comm @@ -333,6 +341,52 @@ def __tosparse_array(self, out=None, is_sparse=False): # Reshape 'allends' to have 9 columns and 'size' rows allends = allends.reshape((size, len(endsarr))) + + # Before we begin computing the dot products we need to know which entries of the output vector tmp2 belong to our rank. + if BoS2 == "s": + # We get the start and endpoint for each sublist in tmp2 + starts2 = tmp2.starts + ends2 = tmp2.ends + # We get the dimensions of the StencilVector + npts2 = self.codomain.npts + # We get the number of space we have + nsp2 = 1 + # We get the number of dimensions the StencilVectorSpace has. + ndim2 = self.codomain.ndim + #We build our ranges of iteration + itterables2 = [] + for ii in range(ndim2): + itterables2.append( + range(starts2[ii], ends2[ii]+1)) + elif BoS2 == "b": + # we collect all starts and ends in two big lists + starts2 = [vi.starts for vi in tmp2] + ends2 = [vi.ends for vi in tmp2] + # We collect the dimension of the BlockVector + npts2 = [sp.npts for sp in self.codomain.spaces] + # We get the number of space we have + nsp2 = len(self.codomain.spaces) + # We get the number of dimensions each space has. + ndim2 = [sp.ndim for sp in self.codomain.spaces] + #We build the range of iteration + itterables2 = [] + # since the size of npts changes denpending on h we need to compute a starting point for + # our row index + spoint2 = 0 + spoint2list = [spoint2] + for hh in range(nsp2): + itterables2aux = [] + for ii in range(ndim2[hh]): + itterables2aux.append( + range(starts2[hh][ii], ends2[hh][ii]+1)) + itterables2.append(itterables2aux) + cummulative2 = 1 + for ii in range(ndim2[hh]): + cummulative2 *= npts2[hh][ii] + spoint2 += cummulative2 + spoint2list.append(spoint2) + + currentrank = 0 # Each rank will take care of setting to 1 each one of its entries while all other entries remain zero. while (currentrank < size): @@ -365,16 +419,41 @@ def __tosparse_array(self, out=None, is_sparse=False): # Compute to which column this iteration belongs col = spoint col += np.ravel_multi_index(i, npts[h]) - if is_sparse == False: - out[:, col] = tmp2.toarray() - else: - aux = tmp2.toarray() - # We now need to now which entries on tmp2 are non-zero and store then in our data list - for l in np.where(aux != 0)[0]: - data.append(aux[l]) - colarr.append(col) - row.append(l) + + # Case in which tmp2 is a StencilVector + if BoS2 == "s": + if is_sparse == False: + #We iterate over the entries of tmp2 that belong to our rank + for ii in itertools.product(*itterables2): + if (tmp2[ii] != 0): + out[np.ravel_multi_index( + ii, npts2), col] = tmp2[ii] + else: + #We iterate over the entries of tmp2 that belong to our rank + for ii in itertools.product(*itterables2): + if (tmp2[ii] != 0): + data.append(tmp2[ii]) + colarr.append(col) + row.append( + np.ravel_multi_index(ii, npts2)) + elif BoS2 =="b": + # We iterate over the stencil vectors inside the BlockVector + for hh in range(nsp2): + itterables2aux = itterables2[hh] + if is_sparse == False: + # We iterate over all the tmp2 entries that belong to rank number currentrank + for ii in itertools.product(*itterables2aux): + if (tmp2[hh][ii] != 0): + out[spoint2list[hh]+np.ravel_multi_index( + ii, npts2[hh]), col] = tmp2[hh][ii] + else: + for ii in itertools.product(*itterables2aux): + if (tmp2[hh][ii] != 0): + data.append(tmp2[hh][ii]) + colarr.append(col) + row.append( + spoint2list[hh]+np.ravel_multi_index(ii, npts2[hh])) ################################# if BoS == "b": if (rank == currentrank): @@ -396,8 +475,8 @@ def __tosparse_array(self, out=None, is_sparse=False): return out else: return sparse.csr_matrix((data, (row, colarr)), shape=(numrows, numcols)) - - + + # Function that returns the local matrix corresponding to the linear operator. Returns a scipy.sparse.csr.csr_matrix. def tosparse(self): """ From db79c1767f40a9cabd5536c70f1b98cc9c95dfe9 Mon Sep 17 00:00:00 2001 From: mateomarin97 Date: Tue, 16 Apr 2024 15:19:48 +0000 Subject: [PATCH 06/18] Added pyccel kernels for construction of dense matrices. --- psydac/linalg/basic.py | 115 ++++++++++++++++------- psydac/linalg/kernels/toarray_kernels.py | 75 +++++++++++++++ 2 files changed, 158 insertions(+), 32 deletions(-) create mode 100644 psydac/linalg/kernels/toarray_kernels.py diff --git a/psydac/linalg/basic.py b/psydac/linalg/basic.py index 163b4c3d4..f2ab12459 100644 --- a/psydac/linalg/basic.py +++ b/psydac/linalg/basic.py @@ -8,6 +8,7 @@ import numpy as np import itertools from scipy import sparse +from psydac.linalg.kernels.toarray_kernels import write_out_stencil_dense_3D, write_out_stencil_dense_2D, write_out_stencil_dense_1D, write_out_block_dense_3D, write_out_block_dense_2D, write_out_block_dense_1D __all__ = ('VectorSpace', 'Vector', 'LinearOperator', 'ZeroOperator', 'IdentityOperator', 'ScaledLinearOperator', 'SumLinearOperator', 'ComposedLinearOperator', 'PowerLinearOperator', 'InverseLinearOperator', 'LinearSolver') @@ -348,43 +349,73 @@ def __tosparse_array(self, out=None, is_sparse=False): starts2 = tmp2.starts ends2 = tmp2.ends # We get the dimensions of the StencilVector - npts2 = self.codomain.npts + npts2 = np.array(self.codomain.npts) # We get the number of space we have nsp2 = 1 # We get the number of dimensions the StencilVectorSpace has. ndim2 = self.codomain.ndim #We build our ranges of iteration - itterables2 = [] - for ii in range(ndim2): - itterables2.append( - range(starts2[ii], ends2[ii]+1)) + if (is_sparse == False): + itterables2 = [] + for ii in range(ndim2): + itterables2.append([starts2[ii], ends2[ii]+1]) + #itterables2.append(range(starts2[ii], ends2[ii]+1)) + itterables2 = np.array(itterables2) + #We also get the StencilVector's pads + pds = np.array(tmp2.pads) + else: + itterables2 = [] + for ii in range(ndim2): + itterables2.append( + range(starts2[ii], ends2[ii]+1)) + elif BoS2 == "b": # we collect all starts and ends in two big lists starts2 = [vi.starts for vi in tmp2] ends2 = [vi.ends for vi in tmp2] # We collect the dimension of the BlockVector - npts2 = [sp.npts for sp in self.codomain.spaces] + npts2 = np.array([sp.npts for sp in self.codomain.spaces]) # We get the number of space we have nsp2 = len(self.codomain.spaces) # We get the number of dimensions each space has. ndim2 = [sp.ndim for sp in self.codomain.spaces] - #We build the range of iteration - itterables2 = [] - # since the size of npts changes denpending on h we need to compute a starting point for - # our row index - spoint2 = 0 - spoint2list = [spoint2] - for hh in range(nsp2): - itterables2aux = [] - for ii in range(ndim2[hh]): - itterables2aux.append( - range(starts2[hh][ii], ends2[hh][ii]+1)) - itterables2.append(itterables2aux) - cummulative2 = 1 - for ii in range(ndim2[hh]): - cummulative2 *= npts2[hh][ii] - spoint2 += cummulative2 - spoint2list.append(spoint2) + if (is_sparse == False): + #We also get the BlockVector's pads + pds = np.array([vi.pads for vi in tmp2]) + #We build the range of iteration + itterables2 = [] + # since the size of npts changes denpending on h we need to compute a starting point for + # our row index + spoint2 = 0 + spoint2list = [np.int64(spoint2)] + for hh in range(nsp2): + itterables2aux = [] + for ii in range(ndim2[hh]): + itterables2aux.append( + [starts2[hh][ii], ends2[hh][ii]+1]) + itterables2.append(itterables2aux) + cummulative2 = 1 + for ii in range(ndim2[hh]): + cummulative2 *= npts2[hh][ii] + spoint2 += cummulative2 + spoint2list.append(spoint2) + else: + itterables2 = [] + # since the size of npts changes denpending on h we need to compute a starting point for + # our row index + spoint2 = 0 + spoint2list = [spoint2] + for hh in range(nsp2): + itterables2aux = [] + for ii in range(ndim2[hh]): + itterables2aux.append( + range(starts2[hh][ii], ends2[hh][ii]+1)) + itterables2.append(itterables2aux) + cummulative2 = 1 + for ii in range(ndim2[hh]): + cummulative2 *= npts2[hh][ii] + spoint2 += cummulative2 + spoint2list.append(spoint2) currentrank = 0 @@ -424,10 +455,20 @@ def __tosparse_array(self, out=None, is_sparse=False): if BoS2 == "s": if is_sparse == False: #We iterate over the entries of tmp2 that belong to our rank - for ii in itertools.product(*itterables2): - if (tmp2[ii] != 0): - out[np.ravel_multi_index( - ii, npts2), col] = tmp2[ii] + #The pyccel kernels are tantamount to this for loop. + #for ii in itertools.product(*itterables2): + #if (tmp2[ii] != 0): + #out[np.ravel_multi_index( + #ii, npts2), col] = tmp2[ii] + if (ndim2 == 3): + write_out_stencil_dense_3D(itterables2, tmp2._data, out, npts2, col, pds) + elif (ndim2 == 2): + write_out_stencil_dense_2D(itterables2, tmp2._data, out, npts2, col, pds) + elif (ndim2 == 1): + write_out_stencil_dense_1D(itterables2, tmp2._data, out, npts2, col, pds) + else: + raise Exception("The codomain dimension must be 3, 2 or 1.") + else: #We iterate over the entries of tmp2 that belong to our rank for ii in itertools.product(*itterables2): @@ -439,15 +480,25 @@ def __tosparse_array(self, out=None, is_sparse=False): elif BoS2 =="b": # We iterate over the stencil vectors inside the BlockVector for hh in range(nsp2): - itterables2aux = itterables2[hh] + if is_sparse == False: + itterables2aux = np.array(itterables2[hh]) # We iterate over all the tmp2 entries that belong to rank number currentrank - for ii in itertools.product(*itterables2aux): - if (tmp2[hh][ii] != 0): - out[spoint2list[hh]+np.ravel_multi_index( - ii, npts2[hh]), col] = tmp2[hh][ii] + #for ii in itertools.product(*itterables2aux): + #if (tmp2[hh][ii] != 0): + #out[spoint2list[hh]+np.ravel_multi_index( + #ii, npts2[hh]), col] = tmp2[hh][ii] + if (ndim2[hh] == 3): + write_out_block_dense_3D(itterables2aux, tmp2[hh]._data, out, npts2[hh], col, pds[hh], spoint2list[hh]) + elif (ndim2[hh] == 2): + write_out_block_dense_2D(itterables2aux, tmp2[hh]._data, out, npts2[hh], col, pds[hh], spoint2list[hh]) + elif (ndim2[hh] == 1): + write_out_block_dense_1D(itterables2aux, tmp2[hh]._data, out, npts2[hh], col, pds[hh], spoint2list[hh]) + else: + raise Exception("The codomain dimension must be 3, 2 or 1.") else: + itterables2aux = itterables2[hh] for ii in itertools.product(*itterables2aux): if (tmp2[hh][ii] != 0): data.append(tmp2[hh][ii]) diff --git a/psydac/linalg/kernels/toarray_kernels.py b/psydac/linalg/kernels/toarray_kernels.py new file mode 100644 index 000000000..78cc2dda1 --- /dev/null +++ b/psydac/linalg/kernels/toarray_kernels.py @@ -0,0 +1,75 @@ +from pyccel.decorators import template + + +#@template(name='T', types=['float[:]', 'complex[:]']) +#@template(name='Tarray', types=['float[:,:]', 'complex[:,:]']) +def write_out_stencil_dense_3D(itterables2:'int[:,:]', tmp2:'float[:,:,:]', out:'float[:,:]', npts2:'int[:]', col:'int64', pds:'int[:]'): + counter0 = 0 + for i0 in range(itterables2[0][0],itterables2[0][1]): + counter1 = 0 + for i1 in range(itterables2[1][0],itterables2[1][1]): + counter2 = 0 + for i2 in range(itterables2[2][0],itterables2[2][1]): + if(tmp2[pds[0]+counter0, pds[1]+counter1, pds[2]+counter2] != 0): + row = i2 + i1*npts2[2] + i0*npts2[2]*npts2[1] + out[row,col] = tmp2[pds[0]+counter0, pds[1]+counter1, pds[2]+counter2] + counter2 += 1 + counter1 += 1 + counter0 += 1 + + +def write_out_stencil_dense_2D(itterables2:'int[:,:]', tmp2:'float[:,:]', out:'float[:,:]', npts2:'int[:]', col:'int64', pds:'int[:]'): + counter0 = 0 + for i0 in range(itterables2[0][0],itterables2[0][1]): + counter1 = 0 + for i1 in range(itterables2[1][0],itterables2[1][1]): + if(tmp2[pds[0]+counter0, pds[1]+counter1] != 0): + row = i1 + i0*npts2[1] + out[row,col] = tmp2[pds[0]+counter0, pds[1]+counter1] + counter1 += 1 + counter0 += 1 + + +def write_out_stencil_dense_1D(itterables2:'int[:,:]', tmp2:'float[:]', out:'float[:,:]', npts2:'int[:]', col:'int64', pds:'int[:]'): + counter0 = 0 + for i0 in range(itterables2[0][0],itterables2[0][1]): + if(tmp2[pds[0]+counter0] != 0): + row = i0 + out[row,col] = tmp2[pds[0]+counter0] + counter0 += 1 + + +def write_out_block_dense_3D(itterables2:'int[:,:]', tmp2:'float[:,:,:]', out:'float[:,:]', npts2:'int[:]', col:'int64', pds:'int[:]', spoint2list:'int64'): + counter0 = 0 + for i0 in range(itterables2[0][0],itterables2[0][1]): + counter1 = 0 + for i1 in range(itterables2[1][0],itterables2[1][1]): + counter2 = 0 + for i2 in range(itterables2[2][0],itterables2[2][1]): + if(tmp2[pds[0]+counter0, pds[1]+counter1, pds[2]+counter2] != 0): + row = i2 + i1*npts2[2] + i0*npts2[2]*npts2[1] + out[spoint2list+row,col] = tmp2[pds[0]+counter0, pds[1]+counter1, pds[2]+counter2] + counter2 += 1 + counter1 += 1 + counter0 += 1 + + +def write_out_block_dense_2D(itterables2:'int[:,:]', tmp2:'float[:,:]', out:'float[:,:]', npts2:'int[:]', col:'int64', pds:'int[:]', spoint2list:'int64'): + counter0 = 0 + for i0 in range(itterables2[0][0],itterables2[0][1]): + counter1 = 0 + for i1 in range(itterables2[1][0],itterables2[1][1]): + if(tmp2[pds[0]+counter0, pds[1]+counter1] != 0): + row = i1 + i0*npts2[1] + out[spoint2list+row,col] = tmp2[pds[0]+counter0, pds[1]+counter1] + counter1 += 1 + counter0 += 1 + + +def write_out_block_dense_1D(itterables2:'int[:,:]', tmp2:'float[:]', out:'float[:,:]', npts2:'int[:]', col:'int64', pds:'int[:]', spoint2list:'int64'): + counter0 = 0 + for i0 in range(itterables2[0][0],itterables2[0][1]): + if(tmp2[pds[0]+counter0] != 0): + row = i0 + out[spoint2list+row,col] = tmp2[pds[0]+counter0] + counter0 += 1 \ No newline at end of file From 4f1a2187093374e867ffc1f7c3fd4d31a4f4c96c Mon Sep 17 00:00:00 2001 From: mateomarin97 Date: Sat, 28 Dec 2024 15:41:11 +0000 Subject: [PATCH 07/18] Parallel tests --- psydac/linalg/tests/toarray_parallel.py | 239 ++++++++++++++++++++++++ 1 file changed, 239 insertions(+) create mode 100644 psydac/linalg/tests/toarray_parallel.py diff --git a/psydac/linalg/tests/toarray_parallel.py b/psydac/linalg/tests/toarray_parallel.py new file mode 100644 index 000000000..1d79d788a --- /dev/null +++ b/psydac/linalg/tests/toarray_parallel.py @@ -0,0 +1,239 @@ +import numpy as np +import scipy.sparse as spa + +from psydac.linalg.direct_solvers import SparseSolver +from psydac.linalg.stencil import StencilVectorSpace, StencilVector, StencilMatrix +from psydac.linalg.block import BlockVectorSpace, BlockVector +from psydac.linalg.block import BlockLinearOperator +from psydac.linalg.utilities import array_to_psydac +from psydac.linalg.kron import KroneckerLinearSolver +from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL +from psydac.ddm.cart import DomainDecomposition, CartDecomposition + +#=============================================================================== +def compare_arrays(arr_psy, arr, rank, atol=1e-14, verbose=False): + '''Assert equality of distributed psydac array and corresponding fraction of cloned numpy array. + Arrays can be block-structured as nested lists/tuples. + + Parameters + ---------- + arr_psy : psydac object + Stencil/Block Vector/Matrix. + + arr : array + Numpy array with same tuple/list structure as arr_psy. If arr is a matrix it can be stencil or band format. + + rank : int + Rank of mpi process + + atol : float + Absolute tolerance used in numpy.allclose() + ''' + if isinstance(arr_psy, StencilVector): + + s = arr_psy.starts + e = arr_psy.ends + + tmp1 = arr_psy[s[0]: e[0] + 1, s[1]: e[1] + 1] + + if arr.ndim == 2: + tmp2 = arr[s[0]: e[0] + 1, s[1]: e[1] + 1] + else: + tmp2 = arr.reshape(arr_psy.space.npts[0], + arr_psy.space.npts[1])[s[0]: e[0] + 1, s[1]: e[1] + 1] + + assert np.allclose(tmp1, tmp2, atol=atol) + + elif isinstance(arr_psy, BlockVector): + + if not (isinstance(arr, tuple) or isinstance(arr, list)): + arrs = np.split(arr, [arr_psy.blocks[0].shape[0], + arr_psy.blocks[0].shape[0] + arr_psy.blocks[1].shape[0]]) + else: + arrs = arr + + for vec_psy, vec in zip(arr_psy.blocks, arrs): + s = vec_psy.starts + e = vec_psy.ends + + tmp1 = vec_psy[s[0]: e[0] + 1, s[1]: e[1] + 1] + + if vec.ndim == 2: + tmp2 = vec[s[0]: e[0] + 1, s[1]: e[1] + 1] + else: + tmp2 = vec.reshape(vec_psy.space.npts[0], vec_psy.space.npts[1])[ + s[0]: e[0] + 1, s[1]: e[1] + 1] + + assert np.allclose(tmp1, tmp2, atol=atol) + else: + raise AssertionError('Wrong input type.') + + if verbose: + print( + f'Rank {rank}: Assertion for array comparison passed with atol={atol}.') + + +def compute_global_starts_ends(domain_decomposition, npts): + global_starts = [None]*len(npts) + global_ends = [None]*len(npts) + + for axis in range(len(npts)): + es = domain_decomposition.global_element_starts[axis] + ee = domain_decomposition.global_element_ends [axis] + + global_ends [axis] = ee.copy() + global_ends [axis][-1] = npts[axis]-1 + global_starts[axis] = np.array([0] + (global_ends[axis][:-1]+1).tolist()) + + return global_starts, global_ends + +def create_equal_random_arrays(W, seedv =123): + np.random.seed(seedv) + arr = [] + arr_psy = BlockVector(W) + + for d, block in enumerate(arr_psy.blocks): + + dims = W.spaces[d].npts + + arr += [np.random.rand(*dims)] + + s = block.starts + e = block.ends + + arr_psy[d][s[0]:e[0] + 1, s[1]:e[1] + 1] = arr[-1][s[0]:e[0] + 1, s[1]:e[1] + 1] + arr_psy.update_ghost_regions() + + return arr, arr_psy + + + +dtype = "float" +n1 = 3 +n2 = 3 +p1 = 1 +p2 = 1 +P1 = True +P2 = True + +# Define a factor for the data +if dtype == "complex": + factor = 1j +else: + factor = 1 + +# set seed for reproducibility +np.random.seed(n1*n2*p1*p2) + +from mpi4py import MPI + +comm = MPI.COMM_WORLD +rank = comm.Get_rank() +size = comm.Get_size() + +D = DomainDecomposition([n1,n2], periods=[P1,P2], comm=comm) + +# Partition the points +npts = [n1,n2] +global_starts, global_ends = compute_global_starts_ends(D, npts) + +cart = CartDecomposition(D, npts, global_starts, global_ends, pads=[p1,p2], shifts=[1,1]) + +# Create vector space, stencil matrix, and stencil vector +V = StencilVectorSpace( cart, dtype=dtype ) +M1 = StencilMatrix( V, V ) +M2 = StencilMatrix( V, V ) +M3 = StencilMatrix( V, V ) +M4 = StencilMatrix( V, V ) +x1 = StencilVector( V ) +x2 = StencilVector( V ) + +s1,s2 = V.starts +e1,e2 = V.ends + +# Fill in stencil matrix values based on diagonal index (periodic!) +for k1 in range(-p1,p1+1): + for k2 in range(-p2,p2+1): + M1[:,:,k1,k2] = factor*k1+k2+10. + M2[:,:,k1,k2] = factor*2.*k1+k2 + M3[:,:,k1,k2] = factor*5*k1+k2 + M4[:,:,k1,k2] = factor*10*k1+k2 + +# If any dimension is not periodic, set corresponding periodic corners to zero +M1.remove_spurious_entries() +M2.remove_spurious_entries() +M3.remove_spurious_entries() +M4.remove_spurious_entries() + +# Fill in vector with random values, then update ghost regions +for i1 in range(s1,e1+1): + for i2 in range(s2,e2+1): + x1[i1,i2] = 2.0*factor*np.random.random() + 1.0 + x2[i1,i2] = 5.0*factor*np.random.random() - 1.0 +x1.update_ghost_regions() +x2.update_ghost_regions() + +# Create and Fill Block objects +W = BlockVectorSpace(V, V) +L = BlockLinearOperator( W, W ) +L[0,0] = M1 +L[0,1] = M2 +L[1,0] = M3 +L[1,1] = M4 + +X = BlockVector(W) +X[0] = x1 +X[1] = x2 + +v1arr1, v1 = create_equal_random_arrays(W, seedv=4568) +v1arr = [] +for i in v1arr1: + aux = i.flatten() + for j in aux: + v1arr.append(j) + +# Test copy with an out +# Create random matrix +N1 = StencilMatrix( V, V ) +N2 = StencilMatrix( V, V ) +N3 = StencilMatrix( V, V ) +N4 = StencilMatrix( V, V ) + +for k1 in range(-p1,p1+1): + for k2 in range(-p2,p2+1): + N1[:,:,k1,k2] = factor*np.random.random() + N2[:,:,k1,k2] = factor*np.random.random() + N3[:,:,k1,k2] = factor*np.random.random() + N4[:,:,k1,k2] = factor*np.random.random() + +K = BlockLinearOperator( W, W ) + +K[0,0] = N1 +K[0,1] = N2 +K[1,0] = N3 +K[1,1] = N4 + +#We compose the two operators +KL = K @ L +KLarr = KL.toarray() + +resultarr = np.matmul(KLarr,v1arr) +resultblock = KL.dot(v1).toarray() +compare_arrays(KL.dot(v1), np.matmul(KLarr, v1arr), rank) +#count = 0 +#while(count < size): + #if(count == rank): + #print("######################",flush=True) + #print("rank = "+str(rank),flush= True) + #print("resultarr : ", flush=True) + #print(resultarr, flush=True) + #print("resultblock :",flush = True) + #print(resultblock, flush = True) + #print("v1arr: " , flush= True) + #print(v1arr, flush= True) + #print("v1.toarray() :", flush= True) + #print(v1.toarray(), flush=True) + #print("######################",flush=True) + #comm.Barrier() + #count += 1 +print("Yesssssss") \ No newline at end of file From b6f5ed255a201a75c4c8528fbd972270f6435320 Mon Sep 17 00:00:00 2001 From: Nathan Date: Mon, 17 Mar 2025 15:26:02 +0100 Subject: [PATCH 08/18] Made the requested clean up of the code. --- psydac/linalg/basic.py | 204 +++++++++++++++++++---------------------- 1 file changed, 94 insertions(+), 110 deletions(-) diff --git a/psydac/linalg/basic.py b/psydac/linalg/basic.py index 59a840d62..d4b8213dd 100644 --- a/psydac/linalg/basic.py +++ b/psydac/linalg/basic.py @@ -280,8 +280,31 @@ def __tosparse_array(self, out=None, is_sparse=False): #if the domain has the cart atrribute, in which case it will be a StencilVectorSpace. if hasattr(self.domain, 'spaces'): BoS = "b" + comm = self.domain.spaces[0].cart.comm + + # we collect all starts and ends in two big lists + starts = [vi.starts for vi in v] + ends = [vi.ends for vi in v] + # We collect the dimension of the BlockVector + npts = [sp.npts for sp in self.domain.spaces] + # We get the number of space we have + nsp = len(self.domain.spaces) + # We get the number of dimensions each space has. + ndim = [sp.ndim for sp in self.domain.spaces] + elif hasattr(self.domain, 'cart'): BoS = "s" + comm = self.domain.cart.comm + + # We get the start and endpoint for each sublist in v + starts = [v.starts] + ends = [v.ends] + # We get the dimensions of the StencilVector + npts = [self.domain.npts] + # We get the number of space we have + nsp = 1 + # We get the number of dimensions the StencilVectorSpace has. + ndim = [self.domain.ndim] else: raise Exception( 'The domain of the LinearOperator must be a BlockVectorSpace or a StencilVectorSpace.') @@ -289,16 +312,80 @@ def __tosparse_array(self, out=None, is_sparse=False): #We also need to know if the codomain is a StencilVectorSpace or a BlockVectorSpace if hasattr(self.codomain, 'spaces'): BoS2 = "b" + + # Before we begin computing the dot products we need to know which entries of the output vector tmp2 belong to our rank. + # we collect all starts and ends in two big lists + starts2 = [vi.starts for vi in tmp2] + ends2 = [vi.ends for vi in tmp2] + # We collect the dimension of the BlockVector + npts2 = np.array([sp.npts for sp in self.codomain.spaces]) + # We get the number of space we have + nsp2 = len(self.codomain.spaces) + # We get the number of dimensions each space has. + ndim2 = [sp.ndim for sp in self.codomain.spaces] + #We build the range of iteration + itterables2 = [] + # since the size of npts changes denpending on h we need to compute a starting point for + # our row index + spoint2 = 0 + if (is_sparse == False): + #We also get the BlockVector's pads + pds = np.array([vi.pads for vi in tmp2]) + spoint2list = [np.int64(spoint2)] + for hh in range(nsp2): + itterables2aux = [] + for ii in range(ndim2[hh]): + itterables2aux.append( + [starts2[hh][ii], ends2[hh][ii]+1]) + itterables2.append(itterables2aux) + cummulative2 = 1 + for ii in range(ndim2[hh]): + cummulative2 *= npts2[hh][ii] + spoint2 += cummulative2 + spoint2list.append(spoint2) + else: + spoint2list = [spoint2] + for hh in range(nsp2): + itterables2aux = [] + for ii in range(ndim2[hh]): + itterables2aux.append( + range(starts2[hh][ii], ends2[hh][ii]+1)) + itterables2.append(itterables2aux) + cummulative2 = 1 + for ii in range(ndim2[hh]): + cummulative2 *= npts2[hh][ii] + spoint2 += cummulative2 + spoint2list.append(spoint2) + elif hasattr(self.codomain, 'cart'): BoS2 = "s" + + # Before we begin computing the dot products we need to know which entries of the output vector tmp2 belong to our rank. + # We get the start and endpoint for each sublist in tmp2 + starts2 = tmp2.starts + ends2 = tmp2.ends + # We get the dimensions of the StencilVector + npts2 = np.array(self.codomain.npts) + # We get the number of space we have + nsp2 = 1 + # We get the number of dimensions the StencilVectorSpace has. + ndim2 = self.codomain.ndim + #We build our ranges of iteration + itterables2 = [] + if (is_sparse == False): + for ii in range(ndim2): + itterables2.append([starts2[ii], ends2[ii]+1]) + itterables2 = np.array(itterables2) + #We also get the StencilVector's pads + pds = np.array(tmp2.pads) + else: + for ii in range(ndim2): + itterables2.append( + range(starts2[ii], ends2[ii]+1)) else: raise Exception( 'The codomain of the LinearOperator must be a BlockVectorSpace or a StencilVectorSpace.') - - if BoS == "b": - comm = self.domain.spaces[0].cart.comm - elif BoS == "s": - comm = self.domain.cart.comm + rank = comm.Get_rank() size = comm.Get_size() @@ -321,29 +408,7 @@ def __tosparse_array(self, out=None, is_sparse=False): data = [] row = [] colarr = [] - - # V is either a BlockVector or a StencilVector depending on the domain of the linear operator. - if BoS == "b": - # we collect all starts and ends in two big lists - starts = [vi.starts for vi in v] - ends = [vi.ends for vi in v] - # We collect the dimension of the BlockVector - npts = [sp.npts for sp in self.domain.spaces] - # We get the number of space we have - nsp = len(self.domain.spaces) - # We get the number of dimensions each space has. - ndim = [sp.ndim for sp in self.domain.spaces] - elif BoS == "s": - # We get the start and endpoint for each sublist in v - starts = [v.starts] - ends = [v.ends] - # We get the dimensions of the StencilVector - npts = [self.domain.npts] - # We get the number of space we have - nsp = 1 - # We get the number of dimensions the StencilVectorSpace has. - ndim = [self.domain.ndim] - + # First each rank is going to need to know the starts and ends of all other ranks startsarr = np.array([starts[i][j] for i in range(nsp) for j in range(ndim[i])], dtype=int) @@ -369,82 +434,6 @@ def __tosparse_array(self, out=None, is_sparse=False): # Reshape 'allends' to have 9 columns and 'size' rows allends = allends.reshape((size, len(endsarr))) - - # Before we begin computing the dot products we need to know which entries of the output vector tmp2 belong to our rank. - if BoS2 == "s": - # We get the start and endpoint for each sublist in tmp2 - starts2 = tmp2.starts - ends2 = tmp2.ends - # We get the dimensions of the StencilVector - npts2 = np.array(self.codomain.npts) - # We get the number of space we have - nsp2 = 1 - # We get the number of dimensions the StencilVectorSpace has. - ndim2 = self.codomain.ndim - #We build our ranges of iteration - if (is_sparse == False): - itterables2 = [] - for ii in range(ndim2): - itterables2.append([starts2[ii], ends2[ii]+1]) - #itterables2.append(range(starts2[ii], ends2[ii]+1)) - itterables2 = np.array(itterables2) - #We also get the StencilVector's pads - pds = np.array(tmp2.pads) - else: - itterables2 = [] - for ii in range(ndim2): - itterables2.append( - range(starts2[ii], ends2[ii]+1)) - - elif BoS2 == "b": - # we collect all starts and ends in two big lists - starts2 = [vi.starts for vi in tmp2] - ends2 = [vi.ends for vi in tmp2] - # We collect the dimension of the BlockVector - npts2 = np.array([sp.npts for sp in self.codomain.spaces]) - # We get the number of space we have - nsp2 = len(self.codomain.spaces) - # We get the number of dimensions each space has. - ndim2 = [sp.ndim for sp in self.codomain.spaces] - if (is_sparse == False): - #We also get the BlockVector's pads - pds = np.array([vi.pads for vi in tmp2]) - #We build the range of iteration - itterables2 = [] - # since the size of npts changes denpending on h we need to compute a starting point for - # our row index - spoint2 = 0 - spoint2list = [np.int64(spoint2)] - for hh in range(nsp2): - itterables2aux = [] - for ii in range(ndim2[hh]): - itterables2aux.append( - [starts2[hh][ii], ends2[hh][ii]+1]) - itterables2.append(itterables2aux) - cummulative2 = 1 - for ii in range(ndim2[hh]): - cummulative2 *= npts2[hh][ii] - spoint2 += cummulative2 - spoint2list.append(spoint2) - else: - itterables2 = [] - # since the size of npts changes denpending on h we need to compute a starting point for - # our row index - spoint2 = 0 - spoint2list = [spoint2] - for hh in range(nsp2): - itterables2aux = [] - for ii in range(ndim2[hh]): - itterables2aux.append( - range(starts2[hh][ii], ends2[hh][ii]+1)) - itterables2.append(itterables2aux) - cummulative2 = 1 - for ii in range(ndim2[hh]): - cummulative2 *= npts2[hh][ii] - spoint2 += cummulative2 - spoint2list.append(spoint2) - - currentrank = 0 # Each rank will take care of setting to 1 each one of its entries while all other entries remain zero. while (currentrank < size): @@ -512,6 +501,7 @@ def __tosparse_array(self, out=None, is_sparse=False): if is_sparse == False: itterables2aux = np.array(itterables2[hh]) # We iterate over all the tmp2 entries that belong to rank number currentrank + #The pyccel kernels are tantamount to this for loop. #for ii in itertools.product(*itterables2aux): #if (tmp2[hh][ii] != 0): #out[spoint2list[hh]+np.ravel_multi_index( @@ -1503,12 +1493,6 @@ def dot(self, v, out=None): self._dot(v).copy(out=out) return out - - def toarray(self): - raise NotImplementedError('toarray() is not defined for MatrixFreeLinearOperator.') - - def tosparse(self): - raise NotImplementedError('tosparse() is not defined for MatrixFreeLinearOperator.') def transpose(self, conjugate=False): if self._dot_transpose is None: From 9b8dda16850e5d68707196073799e1f6c047789d Mon Sep 17 00:00:00 2001 From: mateomarin97 Date: Wed, 19 Mar 2025 12:46:48 +0100 Subject: [PATCH 09/18] Delete psydac/linalg/tests/toarray_parallel.py Erased unnecessary file. --- psydac/linalg/tests/toarray_parallel.py | 239 ------------------------ 1 file changed, 239 deletions(-) delete mode 100644 psydac/linalg/tests/toarray_parallel.py diff --git a/psydac/linalg/tests/toarray_parallel.py b/psydac/linalg/tests/toarray_parallel.py deleted file mode 100644 index 1d79d788a..000000000 --- a/psydac/linalg/tests/toarray_parallel.py +++ /dev/null @@ -1,239 +0,0 @@ -import numpy as np -import scipy.sparse as spa - -from psydac.linalg.direct_solvers import SparseSolver -from psydac.linalg.stencil import StencilVectorSpace, StencilVector, StencilMatrix -from psydac.linalg.block import BlockVectorSpace, BlockVector -from psydac.linalg.block import BlockLinearOperator -from psydac.linalg.utilities import array_to_psydac -from psydac.linalg.kron import KroneckerLinearSolver -from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL -from psydac.ddm.cart import DomainDecomposition, CartDecomposition - -#=============================================================================== -def compare_arrays(arr_psy, arr, rank, atol=1e-14, verbose=False): - '''Assert equality of distributed psydac array and corresponding fraction of cloned numpy array. - Arrays can be block-structured as nested lists/tuples. - - Parameters - ---------- - arr_psy : psydac object - Stencil/Block Vector/Matrix. - - arr : array - Numpy array with same tuple/list structure as arr_psy. If arr is a matrix it can be stencil or band format. - - rank : int - Rank of mpi process - - atol : float - Absolute tolerance used in numpy.allclose() - ''' - if isinstance(arr_psy, StencilVector): - - s = arr_psy.starts - e = arr_psy.ends - - tmp1 = arr_psy[s[0]: e[0] + 1, s[1]: e[1] + 1] - - if arr.ndim == 2: - tmp2 = arr[s[0]: e[0] + 1, s[1]: e[1] + 1] - else: - tmp2 = arr.reshape(arr_psy.space.npts[0], - arr_psy.space.npts[1])[s[0]: e[0] + 1, s[1]: e[1] + 1] - - assert np.allclose(tmp1, tmp2, atol=atol) - - elif isinstance(arr_psy, BlockVector): - - if not (isinstance(arr, tuple) or isinstance(arr, list)): - arrs = np.split(arr, [arr_psy.blocks[0].shape[0], - arr_psy.blocks[0].shape[0] + arr_psy.blocks[1].shape[0]]) - else: - arrs = arr - - for vec_psy, vec in zip(arr_psy.blocks, arrs): - s = vec_psy.starts - e = vec_psy.ends - - tmp1 = vec_psy[s[0]: e[0] + 1, s[1]: e[1] + 1] - - if vec.ndim == 2: - tmp2 = vec[s[0]: e[0] + 1, s[1]: e[1] + 1] - else: - tmp2 = vec.reshape(vec_psy.space.npts[0], vec_psy.space.npts[1])[ - s[0]: e[0] + 1, s[1]: e[1] + 1] - - assert np.allclose(tmp1, tmp2, atol=atol) - else: - raise AssertionError('Wrong input type.') - - if verbose: - print( - f'Rank {rank}: Assertion for array comparison passed with atol={atol}.') - - -def compute_global_starts_ends(domain_decomposition, npts): - global_starts = [None]*len(npts) - global_ends = [None]*len(npts) - - for axis in range(len(npts)): - es = domain_decomposition.global_element_starts[axis] - ee = domain_decomposition.global_element_ends [axis] - - global_ends [axis] = ee.copy() - global_ends [axis][-1] = npts[axis]-1 - global_starts[axis] = np.array([0] + (global_ends[axis][:-1]+1).tolist()) - - return global_starts, global_ends - -def create_equal_random_arrays(W, seedv =123): - np.random.seed(seedv) - arr = [] - arr_psy = BlockVector(W) - - for d, block in enumerate(arr_psy.blocks): - - dims = W.spaces[d].npts - - arr += [np.random.rand(*dims)] - - s = block.starts - e = block.ends - - arr_psy[d][s[0]:e[0] + 1, s[1]:e[1] + 1] = arr[-1][s[0]:e[0] + 1, s[1]:e[1] + 1] - arr_psy.update_ghost_regions() - - return arr, arr_psy - - - -dtype = "float" -n1 = 3 -n2 = 3 -p1 = 1 -p2 = 1 -P1 = True -P2 = True - -# Define a factor for the data -if dtype == "complex": - factor = 1j -else: - factor = 1 - -# set seed for reproducibility -np.random.seed(n1*n2*p1*p2) - -from mpi4py import MPI - -comm = MPI.COMM_WORLD -rank = comm.Get_rank() -size = comm.Get_size() - -D = DomainDecomposition([n1,n2], periods=[P1,P2], comm=comm) - -# Partition the points -npts = [n1,n2] -global_starts, global_ends = compute_global_starts_ends(D, npts) - -cart = CartDecomposition(D, npts, global_starts, global_ends, pads=[p1,p2], shifts=[1,1]) - -# Create vector space, stencil matrix, and stencil vector -V = StencilVectorSpace( cart, dtype=dtype ) -M1 = StencilMatrix( V, V ) -M2 = StencilMatrix( V, V ) -M3 = StencilMatrix( V, V ) -M4 = StencilMatrix( V, V ) -x1 = StencilVector( V ) -x2 = StencilVector( V ) - -s1,s2 = V.starts -e1,e2 = V.ends - -# Fill in stencil matrix values based on diagonal index (periodic!) -for k1 in range(-p1,p1+1): - for k2 in range(-p2,p2+1): - M1[:,:,k1,k2] = factor*k1+k2+10. - M2[:,:,k1,k2] = factor*2.*k1+k2 - M3[:,:,k1,k2] = factor*5*k1+k2 - M4[:,:,k1,k2] = factor*10*k1+k2 - -# If any dimension is not periodic, set corresponding periodic corners to zero -M1.remove_spurious_entries() -M2.remove_spurious_entries() -M3.remove_spurious_entries() -M4.remove_spurious_entries() - -# Fill in vector with random values, then update ghost regions -for i1 in range(s1,e1+1): - for i2 in range(s2,e2+1): - x1[i1,i2] = 2.0*factor*np.random.random() + 1.0 - x2[i1,i2] = 5.0*factor*np.random.random() - 1.0 -x1.update_ghost_regions() -x2.update_ghost_regions() - -# Create and Fill Block objects -W = BlockVectorSpace(V, V) -L = BlockLinearOperator( W, W ) -L[0,0] = M1 -L[0,1] = M2 -L[1,0] = M3 -L[1,1] = M4 - -X = BlockVector(W) -X[0] = x1 -X[1] = x2 - -v1arr1, v1 = create_equal_random_arrays(W, seedv=4568) -v1arr = [] -for i in v1arr1: - aux = i.flatten() - for j in aux: - v1arr.append(j) - -# Test copy with an out -# Create random matrix -N1 = StencilMatrix( V, V ) -N2 = StencilMatrix( V, V ) -N3 = StencilMatrix( V, V ) -N4 = StencilMatrix( V, V ) - -for k1 in range(-p1,p1+1): - for k2 in range(-p2,p2+1): - N1[:,:,k1,k2] = factor*np.random.random() - N2[:,:,k1,k2] = factor*np.random.random() - N3[:,:,k1,k2] = factor*np.random.random() - N4[:,:,k1,k2] = factor*np.random.random() - -K = BlockLinearOperator( W, W ) - -K[0,0] = N1 -K[0,1] = N2 -K[1,0] = N3 -K[1,1] = N4 - -#We compose the two operators -KL = K @ L -KLarr = KL.toarray() - -resultarr = np.matmul(KLarr,v1arr) -resultblock = KL.dot(v1).toarray() -compare_arrays(KL.dot(v1), np.matmul(KLarr, v1arr), rank) -#count = 0 -#while(count < size): - #if(count == rank): - #print("######################",flush=True) - #print("rank = "+str(rank),flush= True) - #print("resultarr : ", flush=True) - #print(resultarr, flush=True) - #print("resultblock :",flush = True) - #print(resultblock, flush = True) - #print("v1arr: " , flush= True) - #print(v1arr, flush= True) - #print("v1.toarray() :", flush= True) - #print(v1.toarray(), flush=True) - #print("######################",flush=True) - #comm.Barrier() - #count += 1 -print("Yesssssss") \ No newline at end of file From cb865492fecd254117988073cdd613c2fa4627f6 Mon Sep 17 00:00:00 2001 From: mateomarin97 Date: Mon, 4 Aug 2025 08:05:14 -0500 Subject: [PATCH 10/18] Create multigrid.py Add multigrid modules --- psydac/linalg/multigrid.py | 1 + 1 file changed, 1 insertion(+) create mode 100644 psydac/linalg/multigrid.py diff --git a/psydac/linalg/multigrid.py b/psydac/linalg/multigrid.py new file mode 100644 index 000000000..8b1378917 --- /dev/null +++ b/psydac/linalg/multigrid.py @@ -0,0 +1 @@ + From cfced7aa9be8a89c7f6853179c468eede49c59c2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Wed, 26 Nov 2025 16:09:59 +0100 Subject: [PATCH 11/18] Use typing.TypeVar in toarray_kernels.py --- psydac/linalg/basic.py | 14 +++++++------- psydac/linalg/kernels/toarray_kernels.py | 17 ++++++++--------- 2 files changed, 15 insertions(+), 16 deletions(-) diff --git a/psydac/linalg/basic.py b/psydac/linalg/basic.py index c90ba6397..198ff1a3b 100644 --- a/psydac/linalg/basic.py +++ b/psydac/linalg/basic.py @@ -285,7 +285,7 @@ def dtype(self): """ - def __tosparse_array(self, out=None, is_sparse=False): + def _tosparse_array(self, out=None, is_sparse=False): """ Transforms the linear operator into a matrix, which is either stored in dense or sparse format. @@ -365,7 +365,7 @@ def __tosparse_array(self, out=None, is_sparse=False): if (is_sparse == False): #We also get the BlockVector's pads pds = np.array([vi.pads for vi in tmp2]) - spoint2list = [np.int64(spoint2)] + spoint2list = [spoint2] for hh in range(nsp2): itterables2aux = [] for ii in range(ndim2[hh]): @@ -374,7 +374,7 @@ def __tosparse_array(self, out=None, is_sparse=False): itterables2.append(itterables2aux) cummulative2 = 1 for ii in range(ndim2[hh]): - cummulative2 *= npts2[hh][ii] + cummulative2 *= npts2[hh][ii].item() spoint2 += cummulative2 spoint2list.append(spoint2) else: @@ -387,7 +387,7 @@ def __tosparse_array(self, out=None, is_sparse=False): itterables2.append(itterables2aux) cummulative2 = 1 for ii in range(ndim2[hh]): - cummulative2 *= npts2[hh][ii] + cummulative2 *= npts2[hh][ii].item() spoint2 += cummulative2 spoint2list.append(spoint2) @@ -499,7 +499,7 @@ def __tosparse_array(self, out=None, is_sparse=False): self.dot(v, out=tmp2) # Compute to which column this iteration belongs col = spoint - col += np.ravel_multi_index(i, npts[h]) + col += int(np.ravel_multi_index(i, npts[h])) # Case in which tmp2 is a StencilVector if BoS2 == "s": @@ -590,7 +590,7 @@ def tosparse(self): The matrix form of the linear operator. If ran in parallel each rank gets the local matrix representation of the linear operator. """ - return self.__tosparse_array(is_sparse=True) + return self._tosparse_array(is_sparse=True) # Function that returns the matrix corresponding to the linear operator. Returns a numpy array. @@ -609,7 +609,7 @@ def toarray(self, out=None): The matrix form of the linear operator. If ran in parallel each rank gets the local matrix representation of the linear operator. """ - return self.__tosparse_array(out=out, is_sparse=False) + return self._tosparse_array(out=out, is_sparse=False) @abstractmethod diff --git a/psydac/linalg/kernels/toarray_kernels.py b/psydac/linalg/kernels/toarray_kernels.py index 78cc2dda1..ee2476880 100644 --- a/psydac/linalg/kernels/toarray_kernels.py +++ b/psydac/linalg/kernels/toarray_kernels.py @@ -1,9 +1,8 @@ -from pyccel.decorators import template +from typing import TypeVar +T = TypeVar('T', float, complex) -#@template(name='T', types=['float[:]', 'complex[:]']) -#@template(name='Tarray', types=['float[:,:]', 'complex[:,:]']) -def write_out_stencil_dense_3D(itterables2:'int[:,:]', tmp2:'float[:,:,:]', out:'float[:,:]', npts2:'int[:]', col:'int64', pds:'int[:]'): +def write_out_stencil_dense_3D(itterables2:'int[:,:]', tmp2:'T[:,:,:]', out:'T[:,:]', npts2:'int[:]', col:int, pds:'int[:]'): counter0 = 0 for i0 in range(itterables2[0][0],itterables2[0][1]): counter1 = 0 @@ -18,7 +17,7 @@ def write_out_stencil_dense_3D(itterables2:'int[:,:]', tmp2:'float[:,:,:]', out: counter0 += 1 -def write_out_stencil_dense_2D(itterables2:'int[:,:]', tmp2:'float[:,:]', out:'float[:,:]', npts2:'int[:]', col:'int64', pds:'int[:]'): +def write_out_stencil_dense_2D(itterables2:'int[:,:]', tmp2:'T[:,:]', out:'T[:,:]', npts2:'int[:]', col:int, pds:'int[:]'): counter0 = 0 for i0 in range(itterables2[0][0],itterables2[0][1]): counter1 = 0 @@ -30,7 +29,7 @@ def write_out_stencil_dense_2D(itterables2:'int[:,:]', tmp2:'float[:,:]', out:'f counter0 += 1 -def write_out_stencil_dense_1D(itterables2:'int[:,:]', tmp2:'float[:]', out:'float[:,:]', npts2:'int[:]', col:'int64', pds:'int[:]'): +def write_out_stencil_dense_1D(itterables2:'int[:,:]', tmp2:'T[:]', out:'T[:,:]', npts2:'int[:]', col:int, pds:'int[:]'): counter0 = 0 for i0 in range(itterables2[0][0],itterables2[0][1]): if(tmp2[pds[0]+counter0] != 0): @@ -39,7 +38,7 @@ def write_out_stencil_dense_1D(itterables2:'int[:,:]', tmp2:'float[:]', out:'flo counter0 += 1 -def write_out_block_dense_3D(itterables2:'int[:,:]', tmp2:'float[:,:,:]', out:'float[:,:]', npts2:'int[:]', col:'int64', pds:'int[:]', spoint2list:'int64'): +def write_out_block_dense_3D(itterables2:'int[:,:]', tmp2:'T[:,:,:]', out:'T[:,:]', npts2:'int[:]', col:int, pds:'int[:]', spoint2list:int): counter0 = 0 for i0 in range(itterables2[0][0],itterables2[0][1]): counter1 = 0 @@ -54,7 +53,7 @@ def write_out_block_dense_3D(itterables2:'int[:,:]', tmp2:'float[:,:,:]', out:'f counter0 += 1 -def write_out_block_dense_2D(itterables2:'int[:,:]', tmp2:'float[:,:]', out:'float[:,:]', npts2:'int[:]', col:'int64', pds:'int[:]', spoint2list:'int64'): +def write_out_block_dense_2D(itterables2:'int[:,:]', tmp2:'T[:,:]', out:'T[:,:]', npts2:'int[:]', col:int, pds:'int[:]', spoint2list:int): counter0 = 0 for i0 in range(itterables2[0][0],itterables2[0][1]): counter1 = 0 @@ -66,7 +65,7 @@ def write_out_block_dense_2D(itterables2:'int[:,:]', tmp2:'float[:,:]', out:'flo counter0 += 1 -def write_out_block_dense_1D(itterables2:'int[:,:]', tmp2:'float[:]', out:'float[:,:]', npts2:'int[:]', col:'int64', pds:'int[:]', spoint2list:'int64'): +def write_out_block_dense_1D(itterables2:'int[:,:]', tmp2:'T[:]', out:'T[:,:]', npts2:'int[:]', col:int, pds:'int[:]', spoint2list:int): counter0 = 0 for i0 in range(itterables2[0][0],itterables2[0][1]): if(tmp2[pds[0]+counter0] != 0): From f2dc22b2668986eff5e6bdd88efbbeea74735ca2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Wed, 26 Nov 2025 16:58:04 +0100 Subject: [PATCH 12/18] Clean up toarray test files --- psydac/linalg/tests/test_toarray.py | 17 +++++---------- psydac/linalg/tests/test_toarray_parallel.py | 23 ++++++-------------- 2 files changed, 13 insertions(+), 27 deletions(-) diff --git a/psydac/linalg/tests/test_toarray.py b/psydac/linalg/tests/test_toarray.py index db29c0257..7068d54ce 100644 --- a/psydac/linalg/tests/test_toarray.py +++ b/psydac/linalg/tests/test_toarray.py @@ -1,21 +1,20 @@ import pytest import numpy as np +from mpi4py import MPI -from psydac.linalg.block import BlockLinearOperator, BlockVector, BlockVectorSpace -from psydac.linalg.basic import LinearOperator, ZeroOperator, IdentityOperator, ComposedLinearOperator, InverseLinearOperator, SumLinearOperator, PowerLinearOperator, ScaledLinearOperator +from psydac.linalg.block import BlockLinearOperator, BlockVector, BlockVectorSpace +from psydac.linalg.basic import InverseLinearOperator, PowerLinearOperator from psydac.linalg.stencil import StencilVectorSpace, StencilVector, StencilMatrix -from psydac.linalg.solvers import ConjugateGradient, inverse +from psydac.linalg.solvers import inverse from psydac.ddm.cart import DomainDecomposition, CartDecomposition -from mpi4py import MPI -from random import random, seed -#=============================================================================== +#=============================================================================== n1array = [2, 7] n2array = [2, 3] p1array = [1, 3] p2array = [1, 3] - +#=============================================================================== def compute_global_starts_ends(domain_decomposition, npts): ndims = len(npts) global_starts = [None]*ndims @@ -97,7 +96,6 @@ def test_square_stencil_basic(n1, n2, p1, p2, P1=False, P2=False): if (P1 or 0 <= i1+k1 < n1) and (P2 or 0 <= i2+k2 < n2): A1[i,j] = nonzero_values[k1,k2] - ### ### Use PowerLinearOperator to test toarray function. ### @@ -127,7 +125,6 @@ def test_square_stencil_basic(n1, n2, p1, p2, P1=False, P2=False): Itest = np.matmul(A1,Sinvarr) assert np.allclose(Itest, np.eye(np.size(Itest[0])),atol= 10**(-5)) - #=============================================================================== @pytest.mark.parametrize('n1', n1array) @@ -184,8 +181,6 @@ def test_square_block_basic(n1, n2, p1, p2, P1=False, P2=False): Bpowarr = Bpow.toarray() assert isinstance(Bpow, PowerLinearOperator) assert np.array_equal(np.matmul(Bpowarr,vbarr),Bpow.dot(vb).toarray()) - - #=============================================================================== # SCRIPT FUNCTIONALITY diff --git a/psydac/linalg/tests/test_toarray_parallel.py b/psydac/linalg/tests/test_toarray_parallel.py index cbe4f77e9..062331b8b 100644 --- a/psydac/linalg/tests/test_toarray_parallel.py +++ b/psydac/linalg/tests/test_toarray_parallel.py @@ -2,17 +2,12 @@ # import pytest import numpy as np -import scipy.sparse as spa - -from psydac.linalg.basic import LinearOperator, ZeroOperator, IdentityOperator, ComposedLinearOperator, InverseLinearOperator, SumLinearOperator, PowerLinearOperator, ScaledLinearOperator -from psydac.linalg.direct_solvers import SparseSolver -from psydac.linalg.stencil import StencilVectorSpace, StencilVector, StencilMatrix -from psydac.linalg.block import BlockVectorSpace, BlockVector -from psydac.linalg.block import BlockLinearOperator -from psydac.linalg.utilities import array_to_psydac -from psydac.linalg.kron import KroneckerLinearSolver -from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL -from psydac.ddm.cart import DomainDecomposition, CartDecomposition + +from psydac.linalg.basic import PowerLinearOperator +from psydac.linalg.stencil import StencilVectorSpace, StencilVector, StencilMatrix +from psydac.linalg.block import BlockVectorSpace, BlockVector +from psydac.linalg.block import BlockLinearOperator +from psydac.ddm.cart import DomainDecomposition, CartDecomposition #=============================================================================== def compare_arrays(arr_psy, arr, rank, atol=1e-14, verbose=False): @@ -125,8 +120,7 @@ def create_equal_random_arrays(W, seedv =123): return arr, arr_psy - - +#=============================================================================== @pytest.mark.parametrize( 'dtype', [float] ) @pytest.mark.parametrize( 'n1', [16, 32] ) @pytest.mark.parametrize( 'n2', [16, 32] ) @@ -227,9 +221,6 @@ def test_block_linear_operator_parallel_dot( dtype, n1, n2, p1, p2, P1, P2 ): assert isinstance(KP, PowerLinearOperator) compare_arrays(KP.dot(v1), np.matmul(KParr, v1arr), rank) - - - #=============================================================================== # SCRIPT FUNCTIONALITY #=============================================================================== From 9d0f4182a1879c0d1ef2c31c7430fb86b8e00b96 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Wed, 26 Nov 2025 17:05:58 +0100 Subject: [PATCH 13/18] Clean up test_tosparse_parallel.py --- psydac/linalg/tests/test_tosparse_parallel.py | 41 ++++++------------- 1 file changed, 13 insertions(+), 28 deletions(-) diff --git a/psydac/linalg/tests/test_tosparse_parallel.py b/psydac/linalg/tests/test_tosparse_parallel.py index 4ea9caf49..36387ae65 100644 --- a/psydac/linalg/tests/test_tosparse_parallel.py +++ b/psydac/linalg/tests/test_tosparse_parallel.py @@ -1,18 +1,11 @@ -# -*- coding: UTF-8 -*- -# import pytest import numpy as np -import scipy.sparse as spa - -from psydac.linalg.basic import LinearOperator, ZeroOperator, IdentityOperator, ComposedLinearOperator, InverseLinearOperator, SumLinearOperator, PowerLinearOperator, ScaledLinearOperator -from psydac.linalg.direct_solvers import SparseSolver -from psydac.linalg.stencil import StencilVectorSpace, StencilVector, StencilMatrix -from psydac.linalg.block import BlockVectorSpace, BlockVector -from psydac.linalg.block import BlockLinearOperator -from psydac.linalg.utilities import array_to_psydac -from psydac.linalg.kron import KroneckerLinearSolver -from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL -from psydac.ddm.cart import DomainDecomposition, CartDecomposition + +from psydac.linalg.basic import PowerLinearOperator +from psydac.linalg.stencil import StencilVectorSpace, StencilVector, StencilMatrix +from psydac.linalg.block import BlockVectorSpace, BlockVector +from psydac.linalg.block import BlockLinearOperator +from psydac.ddm.cart import DomainDecomposition, CartDecomposition #=============================================================================== def compare_arrays(arr_psy, arr, rank, atol=1e-14, verbose=False): @@ -125,8 +118,7 @@ def create_equal_random_arrays(W, seedv =123): return arr, arr_psy - - +#=============================================================================== @pytest.mark.parametrize( 'dtype', [float, complex] ) @pytest.mark.parametrize( 'n1', [16] ) @pytest.mark.parametrize( 'n2', [16, 32] ) @@ -146,7 +138,7 @@ def test_block_linear_operator_parallel_dot( dtype, n1, n2, p1, p2, P1, P2 ): # set seed for reproducibility np.random.seed(n1*n2*p1*p2) - from mpi4py import MPI + from mpi4py import MPI comm = MPI.COMM_WORLD rank = comm.Get_rank() @@ -163,13 +155,12 @@ def test_block_linear_operator_parallel_dot( dtype, n1, n2, p1, p2, P1, P2 ): # Create vector space, stencil matrix, and stencil vector V = StencilVectorSpace( cart, dtype=dtype ) - v0arr , v0 = create_equal_random_arrays(V, seedv=4568) - + v0arr, v0 = create_equal_random_arrays(V, seedv=4568) v0arr = v0arr[0].flatten() - + # Create and Fill Block objects W = BlockVectorSpace(V, V) - + v1arr1, v1 = create_equal_random_arrays(W, seedv=4568) v1arr = [] for i in v1arr1: @@ -196,8 +187,7 @@ def test_block_linear_operator_parallel_dot( dtype, n1, n2, p1, p2, P1, P2 ): K[0,1] = N2 K[1,0] = N3 K[1,1] = N4 - - + ##### #Testing tosparse() for LinearOperators with domain being StencilVectorSpace @@ -227,15 +217,10 @@ def test_block_linear_operator_parallel_dot( dtype, n1, n2, p1, p2, P1, P2 ): assert isinstance(KP, PowerLinearOperator) assert isinstance(KP.domain, BlockVectorSpace) compare_arrays(KP.dot(v1), KLarr.dot(v1arr), rank) - - - - - #=============================================================================== # SCRIPT FUNCTIONALITY #=============================================================================== if __name__ == "__main__": import sys - pytest.main( sys.argv ) \ No newline at end of file + pytest.main(sys.argv) From 47ab655dedc251e1a61bb399ab3523352281f05c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Wed, 26 Nov 2025 17:09:52 +0100 Subject: [PATCH 14/18] Add missing newlines at end of files --- psydac/linalg/kernels/toarray_kernels.py | 2 +- psydac/linalg/tests/test_toarray.py | 2 +- psydac/linalg/tests/test_toarray_parallel.py | 4 +--- 3 files changed, 3 insertions(+), 5 deletions(-) diff --git a/psydac/linalg/kernels/toarray_kernels.py b/psydac/linalg/kernels/toarray_kernels.py index ee2476880..426e2be62 100644 --- a/psydac/linalg/kernels/toarray_kernels.py +++ b/psydac/linalg/kernels/toarray_kernels.py @@ -71,4 +71,4 @@ def write_out_block_dense_1D(itterables2:'int[:,:]', tmp2:'T[:]', out:'T[:,:]', if(tmp2[pds[0]+counter0] != 0): row = i0 out[spoint2list+row,col] = tmp2[pds[0]+counter0] - counter0 += 1 \ No newline at end of file + counter0 += 1 diff --git a/psydac/linalg/tests/test_toarray.py b/psydac/linalg/tests/test_toarray.py index 7068d54ce..9878a10a2 100644 --- a/psydac/linalg/tests/test_toarray.py +++ b/psydac/linalg/tests/test_toarray.py @@ -187,4 +187,4 @@ def test_square_block_basic(n1, n2, p1, p2, P1=False, P2=False): #=============================================================================== if __name__ == "__main__": import sys - pytest.main( sys.argv ) \ No newline at end of file + pytest.main(sys.argv) diff --git a/psydac/linalg/tests/test_toarray_parallel.py b/psydac/linalg/tests/test_toarray_parallel.py index 062331b8b..85639f543 100644 --- a/psydac/linalg/tests/test_toarray_parallel.py +++ b/psydac/linalg/tests/test_toarray_parallel.py @@ -1,5 +1,3 @@ -# -*- coding: UTF-8 -*- -# import pytest import numpy as np @@ -226,4 +224,4 @@ def test_block_linear_operator_parallel_dot( dtype, n1, n2, p1, p2, P1, P2 ): #=============================================================================== if __name__ == "__main__": import sys - pytest.main( sys.argv ) \ No newline at end of file + pytest.main(sys.argv) From 5cb5a72ff3bcaceade808a2925f2f434891662a1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Wed, 26 Nov 2025 17:10:58 +0100 Subject: [PATCH 15/18] Remove unused empty module psydac.linalg.multigrid --- psydac/linalg/multigrid.py | 1 - 1 file changed, 1 deletion(-) delete mode 100644 psydac/linalg/multigrid.py diff --git a/psydac/linalg/multigrid.py b/psydac/linalg/multigrid.py deleted file mode 100644 index 8b1378917..000000000 --- a/psydac/linalg/multigrid.py +++ /dev/null @@ -1 +0,0 @@ - From 2459a739b4f42bc3aed5c81c165e4ccc6bfa27a2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Wed, 26 Nov 2025 17:24:56 +0100 Subject: [PATCH 16/18] Clean up psydac/linalg/basic.py --- psydac/linalg/basic.py | 306 ++++++++++++++++++----------------------- 1 file changed, 136 insertions(+), 170 deletions(-) diff --git a/psydac/linalg/basic.py b/psydac/linalg/basic.py index 198ff1a3b..b6d64ac49 100644 --- a/psydac/linalg/basic.py +++ b/psydac/linalg/basic.py @@ -38,7 +38,6 @@ class VectorSpace(ABC): """ Finite-dimensional vector space V with a scalar (inner) product. - """ @property @abstractmethod @@ -46,7 +45,6 @@ def dimension(self): """ The dimension of a vector space V is the cardinality (i.e. the number of vectors) of a basis of V over its base field. - """ @property @@ -54,7 +52,6 @@ def dimension(self): def dtype(self): """ The data type of the field over which the space is built. - """ @abstractmethod @@ -66,7 +63,6 @@ def zeros(self): ------- null : Vector A new vector object with all components equal to zero. - """ @abstractmethod @@ -98,7 +94,6 @@ def inner(self, x, y): float | complex The scalar product of the two vectors. Note that inner(x, x) is a non-negative real number which is zero if and only if x = 0. - """ @abstractmethod @@ -124,7 +119,6 @@ def axpy(self, a, x, y): class Vector(ABC): """ Element of a vector space V. - """ @property def shape(self): @@ -144,7 +138,6 @@ def inner(self, v): ---------- v : Vector Vector belonging to the same space as self. - """ assert isinstance(v, Vector) assert self.space is v.space @@ -256,7 +249,6 @@ class LinearOperator(ABC): """ Abstract base class for all linear operators acting between two vector spaces V (domain) and W (codomain). - """ @property def shape(self): @@ -281,10 +273,142 @@ def codomain(self): def dtype(self): """ The data type of the field over which the space is built. + """ + + @abstractmethod + def dot(self, v, out=None): + """ Apply the LinearOperator self to the Vector v. + + The result is written to the Vector out, if provided. + + Parameters + ---------- + v : Vector + The vector to which the linear operator (self) is applied. It must + belong to the domain of self. + + out : Vector + The vector in which the result of the operation is stored. It must + belong to the codomain of self. If out is None, a new vector is + created and returned. + Returns + ------- + Vector + The result of the operation. If out is None, a new vector is + returned. Otherwise, the result is stored in out and out is + returned. """ - + @abstractmethod + def transpose(self, conjugate=False): + """ + Transpose the LinearOperator . + + If conjugate is True, return the Hermitian transpose. + """ + + # TODO: check if we should add a copy method!!! + + #------------------------------------- + # Magic methods + #------------------------------------- + def __neg__(self): + """ + Scales itself by -1 and thus returns the addititive inverse as + a new object of the class ScaledLinearOperator. + """ + return ScaledLinearOperator(self.domain, self.codomain, -1.0, self) + + def __mul__(self, c): + """ + Scales a linear operator by a real scalar c by creating an object of the class ScaledLinearOperator, + unless c = 0 or c = 1, in which case either a ZeroOperator or self is returned. + + """ + assert np.isscalar(c) + if c==0: + return ZeroOperator(self.domain, self.codomain) + elif c == 1: + return self + else: + return ScaledLinearOperator(self.domain, self.codomain, c, self) + + def __rmul__(self, c): + """ Calls __mul__ instead. """ + return self * c + + def __matmul__(self, B): + """ + Matrix multiplication using the @ operator. + + If B is a LinearOperator, create a ComposedLinearOperator object. + This is simplified to self if B is an IdentityOperator, and to a + ZeroOperator if B is a ZeroOperator. + + If B is a Vector, the @ operator is treated as a matrix-vector + multiplication and returns the result of self.dot(B). + + Parameters + ---------- + B : LinearOperator | Vector + The object to be multiplied with self. If B is a LinearOperator, + its codomain must be equal to the domain of self. If B is a Vector, + it must belong to the domain of self. + + Returns + ------- + LinearOperator | Vector + If B is a LinearOperator, return a ComposedLinearOperator object, + or a simplification to self or a ZeroOperator. In all cases the + resulting LinearOperator has the same domain as self and the same + codomain as B. If B is a Vector, return the result of self.dot(B), + which is a Vector belonging to the codomain of self. + """ + assert isinstance(B, (LinearOperator, Vector)) + if isinstance(B, LinearOperator): + assert self.domain == B.codomain + if isinstance(B, ZeroOperator): + return ZeroOperator(B.domain, self.codomain) + elif isinstance(B, IdentityOperator): + return self + else: + return ComposedLinearOperator(B.domain, self.codomain, self, B) + else: + return self.dot(B) + + def __add__(self, B): + """ Creates an object of the class SumLinearOperator unless B is a ZeroOperator in which case self is returned. """ + assert isinstance(B, LinearOperator) + if isinstance(B, ZeroOperator): + return self + else: + return SumLinearOperator(self.domain, self.codomain, self, B) + + def __sub__(self, B): + """ Creates an object of the class SumLinearOperator unless B is a ZeroOperator in which case self is returned. """ + assert isinstance(B, LinearOperator) + if isinstance(B, ZeroOperator): + return self + else: + return SumLinearOperator(self.domain, self.codomain, self, -B) + + def __pow__(self, n): + """ Creates an object of class :ref:`PowerLinearOperator `. """ + return PowerLinearOperator(self.domain, self.codomain, self, n) + + def __truediv__(self, c): + """ Divide by scalar. """ + return self * (1.0 / c) + + def __itruediv__(self, c): + """ Divide by scalar, in place. """ + self *= 1.0 / c + return self + + #------------------------------------- + # Methods with default implementation + #------------------------------------- def _tosparse_array(self, out=None, is_sparse=False): """ Transforms the linear operator into a matrix, which is either stored in dense or sparse format. @@ -530,8 +654,6 @@ def _tosparse_array(self, out=None, is_sparse=False): elif BoS2 =="b": # We iterate over the stencil vectors inside the BlockVector for hh in range(nsp2): - - if is_sparse == False: itterables2aux = np.array(itterables2[hh]) # We iterate over all the tmp2 entries that belong to rank number currentrank @@ -572,13 +694,12 @@ def _tosparse_array(self, out=None, is_sparse=False): spoint += cummulative npredim += ndim[h] currentrank += 1 - + if is_sparse == False: return out else: return sparse.csr_matrix((data, (row, colarr)), shape=(numrows, numcols)) - - + # Function that returns the local matrix corresponding to the linear operator. Returns a scipy.sparse.csr.csr_matrix. def tosparse(self): """ @@ -591,8 +712,7 @@ def tosparse(self): matrix representation of the linear operator. """ return self._tosparse_array(is_sparse=True) - - + # Function that returns the matrix corresponding to the linear operator. Returns a numpy array. def toarray(self, out=None): """ @@ -611,142 +731,6 @@ def toarray(self, out=None): """ return self._tosparse_array(out=out, is_sparse=False) - - @abstractmethod - def dot(self, v, out=None): - """ Apply the LinearOperator self to the Vector v. - - The result is written to the Vector out, if provided. - - Parameters - ---------- - v : Vector - The vector to which the linear operator (self) is applied. It must - belong to the domain of self. - - out : Vector - The vector in which the result of the operation is stored. It must - belong to the codomain of self. If out is None, a new vector is - created and returned. - - Returns - ------- - Vector - The result of the operation. If out is None, a new vector is - returned. Otherwise, the result is stored in out and out is - returned. - """ - - @abstractmethod - def transpose(self, conjugate=False): - """ - Transpose the LinearOperator . - - If conjugate is True, return the Hermitian transpose. - """ - - # TODO: check if we should add a copy method!!! - - #------------------------------------- - # Magic methods - #------------------------------------- - def __neg__(self): - """ - Scales itself by -1 and thus returns the addititive inverse as - a new object of the class ScaledLinearOperator. - - """ - return ScaledLinearOperator(self.domain, self.codomain, -1.0, self) - - def __mul__(self, c): - """ - Scales a linear operator by a real scalar c by creating an object of the class ScaledLinearOperator, - unless c = 0 or c = 1, in which case either a ZeroOperator or self is returned. - - """ - assert np.isscalar(c) - if c==0: - return ZeroOperator(self.domain, self.codomain) - elif c == 1: - return self - else: - return ScaledLinearOperator(self.domain, self.codomain, c, self) - - def __rmul__(self, c): - """ Calls __mul__ instead. """ - return self * c - - def __matmul__(self, B): - """ - Matrix multiplication using the @ operator. - - If B is a LinearOperator, create a ComposedLinearOperator object. - This is simplified to self if B is an IdentityOperator, and to a - ZeroOperator if B is a ZeroOperator. - - If B is a Vector, the @ operator is treated as a matrix-vector - multiplication and returns the result of self.dot(B). - - Parameters - ---------- - B : LinearOperator | Vector - The object to be multiplied with self. If B is a LinearOperator, - its codomain must be equal to the domain of self. If B is a Vector, - it must belong to the domain of self. - - Returns - ------- - LinearOperator | Vector - If B is a LinearOperator, return a ComposedLinearOperator object, - or a simplification to self or a ZeroOperator. In all cases the - resulting LinearOperator has the same domain as self and the same - codomain as B. If B is a Vector, return the result of self.dot(B), - which is a Vector belonging to the codomain of self. - """ - assert isinstance(B, (LinearOperator, Vector)) - if isinstance(B, LinearOperator): - assert self.domain == B.codomain - if isinstance(B, ZeroOperator): - return ZeroOperator(B.domain, self.codomain) - elif isinstance(B, IdentityOperator): - return self - else: - return ComposedLinearOperator(B.domain, self.codomain, self, B) - else: - return self.dot(B) - - def __add__(self, B): - """ Creates an object of the class SumLinearOperator unless B is a ZeroOperator in which case self is returned. """ - assert isinstance(B, LinearOperator) - if isinstance(B, ZeroOperator): - return self - else: - return SumLinearOperator(self.domain, self.codomain, self, B) - - def __sub__(self, B): - """ Creates an object of the class SumLinearOperator unless B is a ZeroOperator in which case self is returned. """ - assert isinstance(B, LinearOperator) - if isinstance(B, ZeroOperator): - return self - else: - return SumLinearOperator(self.domain, self.codomain, self, -B) - - def __pow__(self, n): - """ Creates an object of class :ref:`PowerLinearOperator `. """ - return PowerLinearOperator(self.domain, self.codomain, self, n) - - def __truediv__(self, c): - """ Divide by scalar. """ - return self * (1.0 / c) - - def __itruediv__(self, c): - """ Divide by scalar, in place. """ - self *= 1.0 / c - return self - - #------------------------------------- - # Methods with default implementation - #------------------------------------- @property def T(self): """ Calls transpose method to return the transpose of self. """ @@ -829,9 +813,7 @@ def dot_inner(self, v, w): class ZeroOperator(LinearOperator): """ Zero operator mapping any vector from its domain V to the zero vector of its codomain W. - """ - def __new__(cls, domain, codomain=None): assert isinstance(domain, VectorSpace) @@ -928,9 +910,7 @@ class IdentityOperator(LinearOperator): """ Identity operator acting between a vector space V and itself. Useful for example in custom linear operator classes together with the apply_essential_bc method to create projection operators. - """ - def __init__(self, domain, codomain=None): assert isinstance(domain, VectorSpace) @@ -992,9 +972,7 @@ def __matmul__(self, B): class ScaledLinearOperator(LinearOperator): """ A linear operator $A$ scalar multiplied by a constant $c$. - """ - def __init__(self, domain, codomain, c, A): assert isinstance(domain, VectorSpace) @@ -1073,7 +1051,6 @@ def dot(self, v, out=None): class SumLinearOperator(LinearOperator): r""" Sum $\sum_{i=1}^n A_i$ of linear operators $A_1,\dots,A_n$ acting between the same vector spaces V (domain) and W (codomain). - """ def __new__(cls, domain, codomain, *args): @@ -1194,9 +1171,7 @@ def simplify(addends): class ComposedLinearOperator(LinearOperator): r""" Composition $A_n\circ\dots\circ A_1$ of two or more linear operators $A_1,\dots,A_n$. - """ - def __init__(self, domain, codomain, *args): assert isinstance(domain, VectorSpace) @@ -1255,7 +1230,6 @@ def multiplicants(self): r""" A tuple $(A_1,\dots,A_n)$ containing the multiplicants of the linear operator $self = A_n\circ\dots\circ A_1$. - """ return self._multiplicants @@ -1296,7 +1270,6 @@ def dot(self, v, out=None): A = self._multiplicants[0] if out is not None: - A.dot(x, out=out) else: out = A.dot(x) @@ -1314,9 +1287,7 @@ def set_backend(self, backend): class PowerLinearOperator(LinearOperator): r""" Power $A^n$ of a linear operator $A$ for some integer $n\geq 0$. - """ - def __new__(cls, domain, codomain, A, n): assert isinstance(n, int) @@ -1415,7 +1386,6 @@ class InverseLinearOperator(LinearOperator): verbose : bool If True, L2-norm of residual r is printed at each iteration. """ - def __init__(self, A, **kwargs): assert isinstance(A, LinearOperator) @@ -1502,7 +1472,6 @@ def get_options(self, key=None): If `key` is given, get the specific option of interest. If there is no such option, `None` is returned instead. If `key` is not given, get a copy of all the solver options in a dictionary. - """ if key is None: return self._options.copy() @@ -1524,7 +1493,6 @@ def transpose(self, conjugate=False): class LinearSolver(ABC): """ Solver for the square linear system Ax=b, where x and b belong to the same vector space V. - """ @property def shape(self): @@ -1582,9 +1550,7 @@ class MatrixFreeLinearOperator(LinearOperator): # example 2: a truly matrix-free linear operator A = MatrixFreeLinearOperator(domain=V, codomain=V, dot=lambda v: 2*v, dot_transpose=lambda v: 2*v) - """ - def __init__(self, domain, codomain, dot, dot_transpose=None): assert isinstance(domain, VectorSpace) From 3766611e5f3184eaaf6b70a338937e4a1986369e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Wed, 26 Nov 2025 17:41:55 +0100 Subject: [PATCH 17/18] Further clean up psydac/linalg/basic.py --- psydac/linalg/basic.py | 51 +++++++++++++++++++----------------------- 1 file changed, 23 insertions(+), 28 deletions(-) diff --git a/psydac/linalg/basic.py b/psydac/linalg/basic.py index b6d64ac49..8ed52f032 100644 --- a/psydac/linalg/basic.py +++ b/psydac/linalg/basic.py @@ -411,25 +411,25 @@ def __itruediv__(self, c): #------------------------------------- def _tosparse_array(self, out=None, is_sparse=False): """ - Transforms the linear operator into a matrix, which is either stored in dense or sparse format. + Transforms the linear operator into a matrix in dense or sparse format. Parameters ---------- - out : Numpy.ndarray, optional + out : numpy.ndarray, optional If given, the output will be written in-place into this array. is_sparse : bool, optional - If set to True the method returns the matrix as a Scipy sparse matrix, if set to false - it returns the full matrix as a Numpy.ndarray + If True the method returns the matrix as a SciPy sparse matrix, if + False it returns the full matrix as a numpy.ndarray. Returns ------- - out : Numpy.ndarray or scipy.sparse.csr.csr_matrix - The matrix form of the linear operator. If ran in parallel each rank gets the full - matrix representation of the linear operator. + out : numpy.ndarray or scipy.sparse.csr.csr_matrix + The matrix form of the linear operator. If ran in parallel each + rank gets the full matrix representation of the linear operator. """ # v will be the unit vector with which we compute Av = ith column of A. v = self.domain.zeros() - # We define a temporal vector + # We define a temporary vector tmp2 = self.codomain.zeros() #We need to determine if we are a blockvector or a stencilvector but we are not able to use @@ -486,7 +486,7 @@ def _tosparse_array(self, out=None, is_sparse=False): # since the size of npts changes denpending on h we need to compute a starting point for # our row index spoint2 = 0 - if (is_sparse == False): + if not is_sparse: #We also get the BlockVector's pads pds = np.array([vi.pads for vi in tmp2]) spoint2list = [spoint2] @@ -547,7 +547,7 @@ def _tosparse_array(self, out=None, is_sparse=False): rank = comm.Get_rank() size = comm.Get_size() - if (is_sparse == False): + if not is_sparse: if out is None: # We declare the matrix form of our linear operator out = np.zeros( @@ -594,7 +594,7 @@ def _tosparse_array(self, out=None, is_sparse=False): currentrank = 0 # Each rank will take care of setting to 1 each one of its entries while all other entries remain zero. - while (currentrank < size): + while currentrank < size: # since the size of npts changes denpending on h we need to compute a starting point for # our column index spoint = 0 @@ -627,7 +627,7 @@ def _tosparse_array(self, out=None, is_sparse=False): # Case in which tmp2 is a StencilVector if BoS2 == "s": - if is_sparse == False: + if not is_sparse: #We iterate over the entries of tmp2 that belong to our rank #The pyccel kernels are tantamount to this for loop. #for ii in itertools.product(*itterables2): @@ -654,7 +654,7 @@ def _tosparse_array(self, out=None, is_sparse=False): elif BoS2 =="b": # We iterate over the stencil vectors inside the BlockVector for hh in range(nsp2): - if is_sparse == False: + if not is_sparse: itterables2aux = np.array(itterables2[hh]) # We iterate over all the tmp2 entries that belong to rank number currentrank #The pyccel kernels are tantamount to this for loop. @@ -695,39 +695,37 @@ def _tosparse_array(self, out=None, is_sparse=False): npredim += ndim[h] currentrank += 1 - if is_sparse == False: + if not is_sparse: return out else: return sparse.csr_matrix((data, (row, colarr)), shape=(numrows, numcols)) - # Function that returns the local matrix corresponding to the linear operator. Returns a scipy.sparse.csr.csr_matrix. def tosparse(self): """ - Transforms the linear operator into a matrix, which is stored in sparse csr format. + Transforms the linear operator into a matrix in sparse CSR format. Returns ------- - out : Numpy.ndarray or scipy.sparse.csr.csr_matrix - The matrix form of the linear operator. If ran in parallel each rank gets the local - matrix representation of the linear operator. + out : numpy.ndarray or scipy.sparse.csr.csr_matrix + The matrix form of the linear operator. If ran in parallel each + rank gets the local matrix representation of the linear operator. """ return self._tosparse_array(is_sparse=True) - # Function that returns the matrix corresponding to the linear operator. Returns a numpy array. def toarray(self, out=None): """ - Transforms the linear operator into a matrix, which is stored in dense format. + Transforms the linear operator into a matrix in dense format. Parameters ---------- - out : Numpy.ndarray, optional + out : numpy.ndarray, optional If given, the output will be written in-place into this array. Returns ------- - out : Numpy.ndarray - The matrix form of the linear operator. If ran in parallel each rank gets the local - matrix representation of the linear operator. + out : numpy.ndarray + The matrix form of the linear operator. If ran in parallel each + rank gets the local matrix representation of the linear operator. """ return self._tosparse_array(out=out, is_sparse=False) @@ -760,7 +758,6 @@ def idot(self, v, out): out : Vector The vector to be incremented by `self @ v`. It must belong to the codomain of `self`. - """ assert isinstance( v, Vector) assert isinstance(out, Vector) @@ -1213,7 +1210,6 @@ def tmp_vectors(self): """ A tuple containing the storage vectors that are repeatedly being used upon calling the `dot` method. This avoids the creation of new vectors at each call of the `dot` method. - """ return self._tmp_vectors @@ -1425,7 +1421,6 @@ def linop(self): no other options are viable, as it breaks the one-to-one map between the original linear operator $A$ (passed to the constructor) and the current `InverseLinearOperator` object $A^{-1}$. Use with extreme care! - """ return self._A From 612ae7c09823da4284f70a5f9b0fe68241c940a5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Wed, 26 Nov 2025 17:54:26 +0100 Subject: [PATCH 18/18] Further clean up psydac/linalg/basic.py --- psydac/linalg/basic.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/psydac/linalg/basic.py b/psydac/linalg/basic.py index 8ed52f032..500575a47 100644 --- a/psydac/linalg/basic.py +++ b/psydac/linalg/basic.py @@ -530,7 +530,7 @@ def _tosparse_array(self, out=None, is_sparse=False): ndim2 = self.codomain.ndim #We build our ranges of iteration itterables2 = [] - if (is_sparse == False): + if not is_sparse: for ii in range(ndim2): itterables2.append([starts2[ii], ends2[ii]+1]) itterables2 = np.array(itterables2)