Wow, Great answer, thanks, I should probably be able to do it like this, I'll try my best !
I actually do want the assembled matrix now. I have coded the matrix-free version already and it works fine, but I want to use multigrid and I figured that in order to be be most efficient at the coarse level it is best to use -ksp_type preonly -pc_type lu, which requires a true matrix (the matrix I am talking about is the product between a matrix 3dof --> 5dof with a matrix 5dof --> 3dof, which ends up square). But I want to retain a memory light implementation so I use matrix-free formulation at all the levels but the last, which is also light by definition, so I should not suffer much in terms of memory. Thanks Timothée 2015-12-06 22:29 GMT+09:00 Dave May <[email protected]>: > > > On 6 December 2015 at 11:19, Timothée Nicolas <[email protected]> > wrote: > >> Hi, >> >> The answer to the first question is yes. I was puzzled about the second >> part of the answer, but I realize I was not clear enough. I don't want to >> just transfer information between 2 DMDAs, I want to perform an operation >> on vectors with a rectangular matrix. >> > > Okay - I misunderstood. > > > >> To be more explicit, the input vector of the operation is a vector with 3 >> components (hence dof=3), and the results of the operator is put in one >> vector equation + 2 scalar equation, which makes it dof = 5. >> >> So, either I cannot fill my matrix (if the number of rows exceeds what is >> allowed by the DMDA), or I can fill my matrix but then the matrix and >> vectors are not compatible, as in this error : >> > > > Do you need to explicitly assemble this rectangular operator? > Expressing this operator via a matrix-free representation would be > relatively straight forward. > > However, since the DMDA's have the same parallel layout you are in good > shape to define the rectangular matrix if you really want the assembled > thing. Your row space would be defined by the DMDA with block size 5 and > the column space would be defined by the DMDA with block size 3. For > example the matrix would be created like in the following way (assuming a > 2D DMDA) > > PetscInt m,n,M,N; > > > DMDAGetInfo(dm3,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL); > nodes = M *N; > DMDAGetCorners(dm3,NULL,NULL,NULL,&m,&n,NULL); > lnodes = m * n; > MatCreate(PETSC_COMM_WORLD,&A); > MatSetSizes(A,lnodes*5,lnodes*3,nodes*5,nodes*3); > > The sparsity pattern is defined by the connectivity of the DMDA points AND > the different block sizes. > To define the non-zero structure, you need to traverse the part of the > dmda defined by the natural indices given by si <= i < si+m, sj <= j < sj+n > and check whether the the 5x3 block associated with each i,j is living > within the diagonal or off-diagonal block of the matrix. For a given i,j > the off-diagonal is defined by any i',j' associated your stencil which is > not in the range si <= i' <= si + m, sj <= j' < sj + n. (e.g. it is a ghost > node). > > Your arrays counting the non-zero entries on the diagonal (nnz[]) and off > diag (nnz_od[]) can be indexed using i + j * m. This corresponds to the > indexing used by the DMDA. Your non-zero counting function would look > something like the pseudo code below > > PetscInt *nnz,*nnz_od; > > PetscMalloc(sizeof(PetscInt)*m*n*5,&nnz); > PetscMalloc(sizeof(PetscInt)*m*n*5,&nnz_od); > > PetscMemzero(nnz,sizeof(PetscInt)*m*n*5); > PetscMemzero(nnz_od,sizeof(PetscInt)*m*n*5); > > for (j=0; j<n; j++) { > for (i=0; i<m; i++) { > PetscInt idx,d,cnti,cntg; > PetscBool is_ghost; > > idx = i + j * m; > > cnti = 0; > cntg = 0; > for (all_points_in_stencil) { /* User logic to define this loop */ > > is_ghost = PETSC_FALSE; > /*User logic goes here to test if stencil point is a ghost point */ > > /* accumulate counts for the stencil */ > if (is_ghost) { cntg++; } /* values living on a neighbour rank */ > else { cnti++; } /* values local to current rank */ > } > > /* assume all dofs in row and col space in block are coupled to each > other */ > for (d=0; d<5; d++) { > nnz[5*idx+d] = 3 * cnti; > nnz_od[5*idx+d] = 3 * cntg; > } > > } > } > > Thanks, > Dave > > >> >> >> [0]PETSC ERROR: --------------------- Error Message >> -------------------------------------------------------------- >> [0]PETSC ERROR: Nonconforming object sizes >> [0]PETSC ERROR: Mat mat,Vec y: global dim 4900 2940 >> [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html >> for trouble shooting. >> [0]PETSC ERROR: Petsc Release Version 3.6.1, Jul, 22, 2015 >> [0]PETSC ERROR: ./miips on a arch-linux2-c-debug named Carl-9000 by >> timothee Sun Dec 6 19:17:20 2015 >> [0]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++ >> --with-fc=gfortran --download-fblaslapack --download-mpich >> [0]PETSC ERROR: #1 MatMult() line 2215 in >> /home/timothee/Documents/petsc-3.6.1/src/mat/interface/matrix.c >> >> So I assume it means the matrix is assumed to be square. Is there a way >> to change this ? >> >> Best >> >> Timothée >> >> >> >> 2015-12-05 20:41 GMT+09:00 Dave May <[email protected]>: >> >>> >>> >>> On 5 December 2015 at 11:15, Timothée Nicolas < >>> [email protected]> wrote: >>> >>>> Hi, >>>> >>>> I need to create a rectangular matrix which takes vector structured by >>>> a DMDA and transforms them to vectors structured with an other DMDA. It >>>> does not look like there is a routine to do this since apparently >>>> DMCreateMatrix assumes square matrix. What is the clean way to do this ? >>>> >>> >>> Do the DMDA's have exactly the same parallel layout and only differ by >>> the number of DOFs? >>> >>> If yes, and you only want to move entries between vectors associated >>> with the two DMDA's, you can perform the transfer between vectors locally >>> (rank-wise) in a very simple manner using the local number of points from >>> DMDA and the known block sizes (DOFs) associated with your DMDA and the >>> size of the subspace you ultimately want to construct. >>> >>> Thanks, >>> Dave >>> >>> >>>> The problem is I have a DMDA with 8 degrees of freedom (dof) and my >>>> matrix sends vectors living in a subspace with 3 dof to the other 5 dof >>>> subspace. I can define a square matrix living in the large, 8 dof subspace, >>>> but this is of course very inefficient in terms of memory and time. >>>> >>>> Thanks >>>> >>>> Best >>>> >>>> Timothée >>>> >>> >>> >> >
