From fe29b5176a2a4ae964d72b5196c31d15ec8fbc6a Mon Sep 17 00:00:00 2001 From: mateomarin97 Date: Wed, 6 Dec 2023 16:48:33 +0000 Subject: [PATCH 1/7] 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 2/7] 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 3/7] 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 4/7] 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 5/7] 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 6/7] 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 7/7] 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