diff --git a/psydac/linalg/basic.py b/psydac/linalg/basic.py index be06e1f02..59a840d62 100644 --- a/psydac/linalg/basic.py +++ b/psydac/linalg/basic.py @@ -12,7 +12,10 @@ from inspect import signature import numpy as np -from scipy.sparse import coo_matrix +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 + from psydac.utilities.utils import is_real @@ -247,16 +250,344 @@ def codomain(self): @abstractmethod def dtype(self): pass + + 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. - @abstractmethod + 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 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'): + BoS = "s" + 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 + 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] + 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))) + + + # Before we begin computing the dot products we need to know which entries of the output vector tmp2 belong to our rank. + if BoS2 == "s": + # We get the start and endpoint for each sublist in tmp2 + starts2 = tmp2.starts + ends2 = tmp2.ends + # We get the dimensions of the StencilVector + npts2 = np.array(self.codomain.npts) + # We get the number of space we have + nsp2 = 1 + # We get the number of dimensions the StencilVectorSpace has. + ndim2 = self.codomain.ndim + #We build our ranges of iteration + if (is_sparse == False): + itterables2 = [] + for ii in range(ndim2): + itterables2.append([starts2[ii], ends2[ii]+1]) + #itterables2.append(range(starts2[ii], ends2[ii]+1)) + itterables2 = np.array(itterables2) + #We also get the StencilVector's pads + pds = np.array(tmp2.pads) + else: + itterables2 = [] + for ii in range(ndim2): + itterables2.append( + range(starts2[ii], ends2[ii]+1)) + + elif BoS2 == "b": + # we collect all starts and ends in two big lists + starts2 = [vi.starts for vi in tmp2] + ends2 = [vi.ends for vi in tmp2] + # We collect the dimension of the BlockVector + npts2 = np.array([sp.npts for sp in self.codomain.spaces]) + # We get the number of space we have + nsp2 = len(self.codomain.spaces) + # We get the number of dimensions each space has. + ndim2 = [sp.ndim for sp in self.codomain.spaces] + if (is_sparse == False): + #We also get the BlockVector's pads + pds = np.array([vi.pads for vi in tmp2]) + #We build the range of iteration + itterables2 = [] + # since the size of npts changes denpending on h we need to compute a starting point for + # our row index + spoint2 = 0 + spoint2list = [np.int64(spoint2)] + for hh in range(nsp2): + itterables2aux = [] + for ii in range(ndim2[hh]): + itterables2aux.append( + [starts2[hh][ii], ends2[hh][ii]+1]) + itterables2.append(itterables2aux) + cummulative2 = 1 + for ii in range(ndim2[hh]): + cummulative2 *= npts2[hh][ii] + spoint2 += cummulative2 + spoint2list.append(spoint2) + else: + itterables2 = [] + # since the size of npts changes denpending on h we need to compute a starting point for + # our row index + spoint2 = 0 + spoint2list = [spoint2] + for hh in range(nsp2): + itterables2aux = [] + for ii in range(ndim2[hh]): + itterables2aux.append( + range(starts2[hh][ii], ends2[hh][ii]+1)) + itterables2.append(itterables2aux) + cummulative2 = 1 + for ii in range(ndim2[hh]): + cummulative2 *= npts2[hh][ii] + spoint2 += cummulative2 + spoint2list.append(spoint2) + + + currentrank = 0 + # Each rank will take care of setting to 1 each one of its entries while all other entries remain zero. + while (currentrank < size): + # 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() + 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]) + + # 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 + #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): + 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): + + + 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] + 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]) + colarr.append(col) + row.append( + spoint2list[hh]+np.ravel_multi_index(ii, npts2[hh])) + ################################# + if BoS == "b": + if (rank == currentrank): + v[h][i] = 0.0 + v[h].update_ghost_regions() + 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 + + 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): - pass + """ + Transforms the linear operator into a matrix, which is stored in sparse csr format. - @abstractmethod - def toarray(self): - """ Convert to Numpy 2D array. """ - pass + 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): + """ + 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. + + Returns + ------- + out : Numpy.ndarray + The matrix form of the linear operator. If ran in parallel each rank gets the local + matrix representation of the linear operator. + """ + return self.__tosparse_array(out=out, is_sparse=False) + @abstractmethod def dot(self, v, out=None): """ Apply linear operator to Vector v. Result is written to Vector out, if provided.""" @@ -800,15 +1131,12 @@ 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] for Mi in mats[1:]: M = M @ Mi - return coo_matrix(M) + return sparse.coo_matrix(M) def transpose(self, conjugate=False): t_multiplicants = () @@ -907,12 +1235,6 @@ def factorial(self): """ Returns the power to which the operator is raised. """ 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) @@ -1030,12 +1352,6 @@ def _check_options(self, **kwargs): elif key == 'verbose': assert isinstance(value, bool), "verbose must be a bool" - 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): """ Returns the previous convergence information. """ return self._info 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 diff --git a/psydac/linalg/tests/test_toarray.py b/psydac/linalg/tests/test_toarray.py new file mode 100644 index 000000000..db29c0257 --- /dev/null +++ b/psydac/linalg/tests/test_toarray.py @@ -0,0 +1,195 @@ +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. 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 isinstance(Spo3, PowerLinearOperator) + 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): + + # 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..cbe4f77e9 --- /dev/null +++ b/psydac/linalg/tests/test_toarray_parallel.py @@ -0,0 +1,238 @@ +# -*- 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 = [] + if isinstance(W, BlockVectorSpace): + 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() + elif isinstance(W, StencilVectorSpace): + arr_psy = StencilVector(W) + + dims = W.npts + + arr += [np.random.rand(*dims)] + + s = arr_psy.starts + e = arr_psy.ends + + 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 + + + +@pytest.mark.parametrize( 'dtype', [float] ) +@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] ) +@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 ) + + v0arr , v0 = create_equal_random_arrays(V, seedv=4568) + + v0arr = v0arr[0].flatten() + + # Create and Fill Block objects + W = BlockVectorSpace(V, V) + + v1arr1, v1 = create_equal_random_arrays(W, seedv=4568) + v1arr = [] + for i in v1arr1: + aux = i.flatten() + for j in aux: + v1arr.append(j) + + # 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 + + ##### + #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), dtype = dtype) + 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 + KParrloc = KP.toarray() + #We compute the global matrix representation of KP + KParr = np.zeros(np.shape(KParrloc), dtype = dtype) + 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..4ea9caf49 --- /dev/null +++ b/psydac/linalg/tests/test_tosparse_parallel.py @@ -0,0 +1,241 @@ +# -*- 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 = [] + if isinstance(W, BlockVectorSpace): + 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() + elif isinstance(W, StencilVectorSpace): + arr_psy = StencilVector(W) + + dims = W.npts + + arr += [np.random.rand(*dims)] + + s = arr_psy.starts + e = arr_psy.ends + + 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 + + + +@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] ) +@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 ) + + v0arr , v0 = create_equal_random_arrays(V, seedv=4568) + + v0arr = v0arr[0].flatten() + + # Create and Fill Block objects + W = BlockVectorSpace(V, V) + + v1arr1, v1 = create_equal_random_arrays(W, seedv=4568) + v1arr = [] + for i in v1arr1: + aux = i.flatten() + for j in aux: + v1arr.append(j) + + # 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 + + + ##### + #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), dtype=dtype) + 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() + #We get the global matrix of the operator + KLarr = np.zeros(np.shape(KLarrloc), dtype=dtype) + 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) + + + + + + +#=============================================================================== +# SCRIPT FUNCTIONALITY +#=============================================================================== +if __name__ == "__main__": + import sys + pytest.main( sys.argv ) \ No newline at end of file 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