[petsc-users] Discretized fields at quadrature points
Hi, I'd like to be able to represent some discretized fields at the quadrature points of a finite element, so I can easily use them with the plex residual/function evaluations. Is it possible to do this with PetscFECreateDefault and some command line options? I think what I need is PETSCSPACEDG for the space, but I'm not sure what type I should take for the dual space? Thanks, Sander
Re: [petsc-users] VecMTDot vs MatMult
They should be pretty similar. So use whatever makes sense for the rest of the code with respect to the "rectangular dense matrix" or collection of vectors. Barry > On Apr 6, 2016, at 12:07 PM, Phanisri Pradeep Pratapa > wrote: > > Hi, > > I am trying to find out the product of a rectangular dense matrix and a > vector. Could you please let me know if forming the rectangular matrix as an > array of Vecs and then using VecMTDot is better (or worse/similar) than doing > MatMult? (In both cases the matrix and the vectors are parallely distributed) > > Thank you, > > Regards, > > Pradeep
[petsc-users] VecMTDot vs MatMult
Hi, I am trying to find out the product of a rectangular dense matrix and a vector. Could you please let me know if forming the rectangular matrix as an array of Vecs and then using VecMTDot is better (or worse/similar) than doing MatMult? (In both cases the matrix and the vectors are parallely distributed) Thank you, Regards, Pradeep
Re: [petsc-users] Prelloc: Get coordinates of local diagonal submatrix
On Wed, Apr 6, 2016 at 10:12 AM, Florian Lindner
wrote:
> Hello,
>
> in order to preallocate I have to know whether a non-zero will be in the
> local diagonal submatrix or not.
>
>
> http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatMPIAIJSetPreallocation.html
>
> says:
>
> The DIAGONAL portion of the local submatrix of a processor can be defined
> as the submatrix which is obtained by extraction the part corresponding to
> the rows r1-r2 and columns c1-c2 of the global matrix, where r1 is the
> first row that belongs to the processor, r2 is the last row belonging to
> the this processor, and c1-c2 is range of indices of the local part of a
> vector suitable for applying the matrix to. This is an mxn matrix. In the
> common case of a square matrix, the row and column ranges are the same and
> the DIAGONAL part is also square. The remaining portion of the local
> submatrix (mxN) constitute the OFF-DIAGONAL portion.
>
> Therefore I compute the begin and end rows and cols like that:
>
> int ownerRangeBeginA = _matrixA.ownerRange().first;
> int ownerRangeEndA = _matrixA.ownerRange().second;
>
> // Uses: MatGetVecs(_matrixA.matrix, &vector, nullptr);
> petsc::Vector diagTest{_matrixA, "test", petsc::Vector::RIGHT};
>
> int localDiagColBegin = diagTest.ownerRange().first;
> int localDiagColEnd = diagTest.ownerRange().second;
>
> Debug("Local Submatrix Rows = " << ownerRangeBeginA << " / " <<
> ownerRangeEndA <<
> ", Local Submatrix Cols = " << localDiagColBegin << " / " <<
> localDiagColEnd);
>
>
> It's a little bit tainted by my PETSc C++ helper function, but I think you
> get the code.
>
> Is there a way to do it more elegantly? Without instantiating a vector
> just for the purpose of gettings its owner range? And then testing a if an
> index is inside the local diagonal submatrix?
>
http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatGetOwnershipRangesColumn.html
Thanks,
Matt
> Thanks,
> Florian
>
--
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
[petsc-users] Prelloc: Get coordinates of local diagonal submatrix
Hello,
in order to preallocate I have to know whether a non-zero will be in the local
diagonal submatrix or not.
http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatMPIAIJSetPreallocation.html
says:
The DIAGONAL portion of the local submatrix of a processor can be defined as
the submatrix which is obtained by extraction the part corresponding to the
rows r1-r2 and columns c1-c2 of the global matrix, where r1 is the first row
that belongs to the processor, r2 is the last row belonging to the this
processor, and c1-c2 is range of indices of the local part of a vector suitable
for applying the matrix to. This is an mxn matrix. In the common case of a
square matrix, the row and column ranges are the same and the DIAGONAL part is
also square. The remaining portion of the local submatrix (mxN) constitute the
OFF-DIAGONAL portion.
Therefore I compute the begin and end rows and cols like that:
int ownerRangeBeginA = _matrixA.ownerRange().first;
int ownerRangeEndA = _matrixA.ownerRange().second;
// Uses: MatGetVecs(_matrixA.matrix, &vector, nullptr);
petsc::Vector diagTest{_matrixA, "test", petsc::Vector::RIGHT};
int localDiagColBegin = diagTest.ownerRange().first;
int localDiagColEnd = diagTest.ownerRange().second;
Debug("Local Submatrix Rows = " << ownerRangeBeginA << " / " <<
ownerRangeEndA <<
", Local Submatrix Cols = " << localDiagColBegin << " / " <<
localDiagColEnd);
It's a little bit tainted by my PETSc C++ helper function, but I think you get
the code.
Is there a way to do it more elegantly? Without instantiating a vector just for
the purpose of gettings its owner range? And then testing a if an index is
inside the local diagonal submatrix?
Thanks,
Florian
Re: [petsc-users] How to extract array portions
Fantastic, it definitely solved my problem! Thank you! Should I always create the vectors with MatCreateVecs when I use them for matrix-vector multiplications? I used to create vectors with VecCreate, and it has always worked fine with "standard" matrices (matrices not created with MATNEST structure). In the example at http://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/examples/tutorials/ex2.c.html I read the following: "PETSc automatically generates appropriately partitioned matrices and vectors when MatCreate() and VecCreate() are used with the same communicator." It is not clear to me when I can simply use VecCreate (as I've always done so far), and when to use MatCreateVecs... On Wed, 6 Apr 2016 14:38:12 +0200 Dave May wrote: On 6 April 2016 at 14:12, FRANCAVILLA MATTEO ALESSANDRO wrote: Thanks Dave, that's exactly what I was looking for, and it (partially) solved my problem. The MatMult operation now works fine, unfortunately not for every number of MPI processes (while the MatMult involving each submatrix works fine for any number of processes): matArray(1) = L_discretized matArray(2) = K_discretized matArray(3) = L_MR matArray(4) = K_MR CALL MatCreateNest(PETSC_COMM_WORLD,2,PETSC_NULL_OBJECT,2,PETSC_NULL_OBJECT,matArray,TheMatrix,IERR) CALL MatSetFromOptions(TheMatrix,IERR) ! CALL VecCreate(PETSC_COMM_WORLD,x,IERR) CALL VecSetSizes(x,PETSC_DECIDE,2*NumberOfUnknowns,IERR) CALL VecSetFromOptions(x,IERR) CALL VecSet(x,0.+0.*PETSC_i,IERR) CALL VecSetValues(x,1,0,1.+0.*PETSC_i,INSERT_VALUES,IERR) CALL VecAssemblyBegin(x,IERR) CALL VecAssemblyEnd(x,IERR) CALL VecCreate(PETSC_COMM_WORLD,y,IERR) CALL VecSetSizes(y,PETSC_DECIDE,NumberOfUnknowns+NumberOfMeasures,IERR) CALL VecSetFromOptions(y,IERR) CALL MatMult(TheMatrix,x,y,IERR) I guess there is something wrong with the distribution of vectors between the MPI processes. With the specific sizes of my problem everything works fine with 1-2-3-6 processes, otherwise I get "Nonconforming object sizes" error, e.g. with 8 processes: [2]PETSC ERROR: - Error Message -- [2]PETSC ERROR: Nonconforming object sizes [2]PETSC ERROR: Mat mat,Vec y: local dim 383 382 [2]PETSC ERROR: Petsc Release Version 3.6.3, Dec, 03, 2015 [2]PETSC ERROR: Nonconforming object sizes ./ISS on a arch-linux2-c-debug named alessandro-HP-ProDesk-490-G2-MT by alessandro Wed Apr 6 13:31:45 2016 [2]PETSC ERROR: Configure options --with-blas-lapack-dir=/opt/intel/composer_xe_2015.2.164/mkl --with-mpi-cc=/home/alessandro/.openmpi-1.8.3/bin/mpicc --with-mpi-f90=/home/alessandro/.openmpi-1.8.3/bin/mpif90 --with-mpiexec=/home/alessandro/.openmpi-1.8.3/bin/mpiexec --with-scalar-type=complex [2]PETSC ERROR: #1 MatMult() line 2216 in /home/alessandro/CODE/PETSc/petsc-3.6.3/src/mat/interface/matrix.c ... [6]PETSC ERROR: Mat mat,Vec y: local dim 381 382 ... [7]PETSC ERROR: Mat mat,Vec y: local dim 381 382 ... [1]PETSC ERROR: Mat mat,Vec y: local dim 383 382 ... I noticed that the code works (runs and provides correct results) when the number of processes divides the length of y2 (number of rows of A21 and A22); however when I apply MatMult separately on A21 and A22 everything works fine with any number of processes Any idea where I'm messing up here? Sounds like the large vector you are using to multiple with the matnest object wasn't created via MatCreateVecs(). If you use this method you will obtain a vector with the correct local and global sizes. Thanks, Alessandro On Wed, 6 Apr 2016 11:16:43 +0200 Dave May wrote: On 6 April 2016 at 11:08, FRANCAVILLA MATTEO ALESSANDRO < d019...@polito.it> wrote: Hi, I'm trying to set a linear problem with a 2x2 block matrix (say A=[A11 A12; A21 A22] in matlab notation). Obviously I could generate a single matrix, but I would like to keep it as a 2x2 block matrix with the 4 blocks generated and used independently (e.g., because at some point I need to change some of those blocks to matrix-free). My plan was then to generate a MatShell and set a MATOP_MULT shell operation to explicitly compute y=A*x by splitting x=[x1;x2] and y=[y1;y2] and carrying out the single matrix-vector multiplications, and reassembling the final result. Take a look at MatNest - it does exactly what you want. http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatCreateNest.html http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MATNEST.html#MATNEST Thanks, Dave My problem is that my efforts in extracting x1 and x2 from the original x vector have failed. The following Fortran code is what seemed to me the most reasonable thing to do to set the matrix-vector multiplication, but it does not work !== subroutine MyMult(A,x,y,IERR) implicit none #include #include #include #include Mat A Vec x,y, x1, x2, y1, y2 I
Re: [petsc-users] How to extract array portions
On 6 April 2016 at 16:20, FRANCAVILLA MATTEO ALESSANDRO wrote: > Fantastic, it definitely solved my problem! Thank you! > > Should I always create the vectors with MatCreateVecs when I use them for > matrix-vector multiplications? I used to create vectors with VecCreate, and > it has always worked fine with "standard" matrices (matrices not created > with MATNEST structure). > Yes > > In the example at > http://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/examples/tutorials/ex2.c.html > I read the following: > > "PETSc automatically generates appropriately partitioned matrices and > vectors when MatCreate() and VecCreate() are used with the same > communicator." > This comment also assumes that you use PETSC_DECIDE for the local sizes. So if you create a matrix with global size M x M and a vector of global size M and you (i) use PETSC_DECIDE for the local size for the row/col with MatCreate() and PETSC_DECIDE for the local row for VecCreate(); and (ii) you use the same communicator - then the row partitioning adopted by the matrix will be consistent with the row partitioning for the vector. > It is not clear to me when I can simply use VecCreate (as I've always done > so far), and when to use MatCreateVecs... > You should always use MatCreateVecs(). The implementation of Mat will give you back a consistent vector. Thanks, Dave > > > > On Wed, 6 Apr 2016 14:38:12 +0200 > > Dave May wrote: > >> On 6 April 2016 at 14:12, FRANCAVILLA MATTEO ALESSANDRO < >> d019...@polito.it> >> wrote: >> >> Thanks Dave, that's exactly what I was looking for, and it (partially) >>> solved my problem. >>> >>> The MatMult operation now works fine, unfortunately not for every number >>> of MPI processes (while the MatMult involving each submatrix works fine >>> for >>> any number of processes): >>> >>> >>> matArray(1) = L_discretized >>> matArray(2) = K_discretized >>> matArray(3) = L_MR >>> matArray(4) = K_MR >>> CALL >>> >>> MatCreateNest(PETSC_COMM_WORLD,2,PETSC_NULL_OBJECT,2,PETSC_NULL_OBJECT,matArray,TheMatrix,IERR) >>> CALL MatSetFromOptions(TheMatrix,IERR) >>> ! >>> CALL VecCreate(PETSC_COMM_WORLD,x,IERR) >>> CALL VecSetSizes(x,PETSC_DECIDE,2*NumberOfUnknowns,IERR) >>> CALL VecSetFromOptions(x,IERR) >>> CALL VecSet(x,0.+0.*PETSC_i,IERR) >>> CALL VecSetValues(x,1,0,1.+0.*PETSC_i,INSERT_VALUES,IERR) >>> CALL VecAssemblyBegin(x,IERR) >>> CALL VecAssemblyEnd(x,IERR) >>> CALL VecCreate(PETSC_COMM_WORLD,y,IERR) >>> CALL VecSetSizes(y,PETSC_DECIDE,NumberOfUnknowns+NumberOfMeasures,IERR) >>> CALL VecSetFromOptions(y,IERR) >>> CALL MatMult(TheMatrix,x,y,IERR) >>> >>> >>> >>> I guess there is something wrong with the distribution of vectors between >>> the MPI processes. With the specific sizes of my problem everything works >>> fine with 1-2-3-6 processes, otherwise I get "Nonconforming object sizes" >>> error, e.g. with 8 processes: >>> >>> [2]PETSC ERROR: - Error Message >>> -- >>> [2]PETSC ERROR: Nonconforming object sizes >>> [2]PETSC ERROR: Mat mat,Vec y: local dim 383 382 >>> [2]PETSC ERROR: Petsc Release Version 3.6.3, Dec, 03, 2015 >>> [2]PETSC ERROR: Nonconforming object sizes >>> ./ISS on a arch-linux2-c-debug named alessandro-HP-ProDesk-490-G2-MT by >>> alessandro Wed Apr 6 13:31:45 2016 >>> [2]PETSC ERROR: Configure options >>> --with-blas-lapack-dir=/opt/intel/composer_xe_2015.2.164/mkl >>> --with-mpi-cc=/home/alessandro/.openmpi-1.8.3/bin/mpicc >>> --with-mpi-f90=/home/alessandro/.openmpi-1.8.3/bin/mpif90 >>> --with-mpiexec=/home/alessandro/.openmpi-1.8.3/bin/mpiexec >>> --with-scalar-type=complex >>> [2]PETSC ERROR: #1 MatMult() line 2216 in >>> /home/alessandro/CODE/PETSc/petsc-3.6.3/src/mat/interface/matrix.c >>> ... >>> [6]PETSC ERROR: Mat mat,Vec y: local dim 381 382 >>> ... >>> [7]PETSC ERROR: Mat mat,Vec y: local dim 381 382 >>> ... >>> [1]PETSC ERROR: Mat mat,Vec y: local dim 383 382 >>> ... >>> >>> >>> I noticed that the code works (runs and provides correct results) when >>> the >>> number of processes divides the length of y2 (number of rows of A21 and >>> A22); however when I apply MatMult separately on A21 and A22 everything >>> works fine with any number of processes >>> >>> Any idea where I'm messing up here? >>> >>> >> Sounds like the large vector you are using to multiple with the matnest >> object wasn't created via MatCreateVecs(). If you use this method you will >> obtain a vector with the correct local and global sizes. >> >> >> >> >> >>> >>> Thanks, >>> Alessandro >>> >>> >>> >>> >>> On Wed, 6 Apr 2016 11:16:43 +0200 >>> Dave May wrote: >>> >>> On 6 April 2016 at 11:08, FRANCAVILLA MATTEO ALESSANDRO < d019...@polito.it> wrote: Hi, > > I'm trying to set a linear problem with a 2x2 block matrix (say A=[A11 > A12; A21 A22] in matlab notation). Obviously I could generate a single > matrix, but I would like to keep it as a 2x2 block matrix with the 4
Re: [petsc-users] SLEPc eigenpairs relative error
> El 6 abr 2016, a las 14:43, Thibault Chevret > escribió: > > Dear all, > > I am interested in computing the eigenpairs of a sparse real non-symmetric > matrix (~8000x8000). I use SLEPc for this purpose, specifying only nev, ncv > and tol as options. > > My problem is, I guess, that I would like to compute the whole spectra, and > I've got some issues regarding the errors on some eigenvalues (ev). > Real ones are OK, their relative errors being below 1e-10, but for some > complex ones I can't get below 1e-1. > > This is greatly improved if I specify -eps_larger_imaginary, but then the > problem is inversed as the errors associated to real ev become worse. > My solution, for now, is to merge results obtained for complex ev using a > resolution with -eps_larger_imaginary, and those obtained for real ev without > this option. This is quite tedious, as I have to resolve my matrix twice, but > it works fine. > > Am I doing something wrong ? I've read on another post that SLEPc should not > be used for the computation of the whole spectra, but is there by any chance > any way to improve the results using only one calculation ? > > Thank you for your time. > > Thibault SLEPc is not intended for computing the whole spectrum. For this you can use LAPACK (or ScaLAPACK) treating the matrix as dense. LAPACK can also provide information about conditioning, which may be causing large residuals. In SLEPc, you can try setting a value for the mpd parameter (e.g. -eps_mpd 400). See explanation in section 2.6.5 of the users manual. But again SLEPc is designed for a subset of the spectrum. Jose
[petsc-users] SLEPc eigenpairs relative error
Dear all, I am interested in computing the eigenpairs of a sparse real non-symmetric matrix (~8000x8000). I use SLEPc for this purpose, specifying only nev, ncv and tol as options. My problem is, I guess, that I would like to compute the whole spectra, and I've got some issues regarding the errors on some eigenvalues (ev). Real ones are OK, their relative errors being below 1e-10, but for some complex ones I can't get below 1e-1. This is greatly improved if I specify -eps_larger_imaginary, but then the problem is inversed as the errors associated to real ev become worse. My solution, for now, is to merge results obtained for complex ev using a resolution with -eps_larger_imaginary, and those obtained for real ev without this option. This is quite tedious, as I have to resolve my matrix twice, but it works fine. Am I doing something wrong ? I've read on another post that SLEPc should not be used for the computation of the whole spectra, but is there by any chance any way to improve the results using only one calculation ? Thank you for your time. Thibault
Re: [petsc-users] How to extract array portions
On 6 April 2016 at 14:12, FRANCAVILLA MATTEO ALESSANDRO wrote: > Thanks Dave, that's exactly what I was looking for, and it (partially) > solved my problem. > > The MatMult operation now works fine, unfortunately not for every number > of MPI processes (while the MatMult involving each submatrix works fine for > any number of processes): > > > matArray(1) = L_discretized > matArray(2) = K_discretized > matArray(3) = L_MR > matArray(4) = K_MR > CALL > MatCreateNest(PETSC_COMM_WORLD,2,PETSC_NULL_OBJECT,2,PETSC_NULL_OBJECT,matArray,TheMatrix,IERR) > CALL MatSetFromOptions(TheMatrix,IERR) > ! > CALL VecCreate(PETSC_COMM_WORLD,x,IERR) > CALL VecSetSizes(x,PETSC_DECIDE,2*NumberOfUnknowns,IERR) > CALL VecSetFromOptions(x,IERR) > CALL VecSet(x,0.+0.*PETSC_i,IERR) > CALL VecSetValues(x,1,0,1.+0.*PETSC_i,INSERT_VALUES,IERR) > CALL VecAssemblyBegin(x,IERR) > CALL VecAssemblyEnd(x,IERR) > CALL VecCreate(PETSC_COMM_WORLD,y,IERR) > CALL VecSetSizes(y,PETSC_DECIDE,NumberOfUnknowns+NumberOfMeasures,IERR) > CALL VecSetFromOptions(y,IERR) > CALL MatMult(TheMatrix,x,y,IERR) > > > > I guess there is something wrong with the distribution of vectors between > the MPI processes. With the specific sizes of my problem everything works > fine with 1-2-3-6 processes, otherwise I get "Nonconforming object sizes" > error, e.g. with 8 processes: > > [2]PETSC ERROR: - Error Message > -- > [2]PETSC ERROR: Nonconforming object sizes > [2]PETSC ERROR: Mat mat,Vec y: local dim 383 382 > [2]PETSC ERROR: Petsc Release Version 3.6.3, Dec, 03, 2015 > [2]PETSC ERROR: Nonconforming object sizes > ./ISS on a arch-linux2-c-debug named alessandro-HP-ProDesk-490-G2-MT by > alessandro Wed Apr 6 13:31:45 2016 > [2]PETSC ERROR: Configure options > --with-blas-lapack-dir=/opt/intel/composer_xe_2015.2.164/mkl > --with-mpi-cc=/home/alessandro/.openmpi-1.8.3/bin/mpicc > --with-mpi-f90=/home/alessandro/.openmpi-1.8.3/bin/mpif90 > --with-mpiexec=/home/alessandro/.openmpi-1.8.3/bin/mpiexec > --with-scalar-type=complex > [2]PETSC ERROR: #1 MatMult() line 2216 in > /home/alessandro/CODE/PETSc/petsc-3.6.3/src/mat/interface/matrix.c > ... > [6]PETSC ERROR: Mat mat,Vec y: local dim 381 382 > ... > [7]PETSC ERROR: Mat mat,Vec y: local dim 381 382 > ... > [1]PETSC ERROR: Mat mat,Vec y: local dim 383 382 > ... > > > I noticed that the code works (runs and provides correct results) when the > number of processes divides the length of y2 (number of rows of A21 and > A22); however when I apply MatMult separately on A21 and A22 everything > works fine with any number of processes > > Any idea where I'm messing up here? > Sounds like the large vector you are using to multiple with the matnest object wasn't created via MatCreateVecs(). If you use this method you will obtain a vector with the correct local and global sizes. > > > Thanks, > Alessandro > > > > > On Wed, 6 Apr 2016 11:16:43 +0200 > Dave May wrote: > >> On 6 April 2016 at 11:08, FRANCAVILLA MATTEO ALESSANDRO < >> d019...@polito.it> >> wrote: >> >> Hi, >>> >>> I'm trying to set a linear problem with a 2x2 block matrix (say A=[A11 >>> A12; A21 A22] in matlab notation). Obviously I could generate a single >>> matrix, but I would like to keep it as a 2x2 block matrix with the 4 >>> blocks >>> generated and used independently (e.g., because at some point I need to >>> change some of those blocks to matrix-free). My plan was then to >>> generate a >>> MatShell and set a MATOP_MULT shell operation to explicitly compute y=A*x >>> by splitting x=[x1;x2] and y=[y1;y2] and carrying out the single >>> matrix-vector multiplications, and reassembling the final result. >>> >> >> >> Take a look at MatNest - it does exactly what you want. >> >> >> http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatCreateNest.html >> >> >> http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MATNEST.html#MATNEST >> >> Thanks, >> Dave >> >> >> My problem is that my efforts in extracting x1 and x2 from the original x >>> vector have failed. The following Fortran code is what seemed to me the >>> most reasonable thing to do to set the matrix-vector multiplication, but >>> it >>> does not work >>> >>> >>> !== >>> subroutine MyMult(A,x,y,IERR) >>> >>> implicit none >>> >>> #include >>> #include >>> #include >>> #include >>> >>> Mat A >>> Vec x,y, x1, x2, y1, y2 >>> IS :: ix1, ix2, iy1, iy2 >>> PetscErrorCode ierr >>> PetscIntM, N, II, NumberOfUnknowns, NumberOfMeasures >>> PetscInt, ALLOCATABLE :: INDICES(:) >>> PetscMPIInt :: RANK >>> >>> CALL VecGetSize(x,N,IERR) >>> NumberOfUnknowns = N / 2 >>> CALL VecGetSize(y,N,IERR) >>> NumberOfMeasures = N - NumberOfUnknowns >>> >>> N = NumberOfUnknowns + MAX(NumberOfUnknowns,NumberOfMeasures) >>> ALLOCATE(INDICES(N)) >>> INDICES = (/ (ii, ii=0, N-1) /) >>> >>> CALL >>> >>> ISCreateGener
Re: [petsc-users] How to extract array portions
Thanks Dave, that's exactly what I was looking for, and it (partially) solved my problem. The MatMult operation now works fine, unfortunately not for every number of MPI processes (while the MatMult involving each submatrix works fine for any number of processes): matArray(1) = L_discretized matArray(2) = K_discretized matArray(3) = L_MR matArray(4) = K_MR CALL MatCreateNest(PETSC_COMM_WORLD,2,PETSC_NULL_OBJECT,2,PETSC_NULL_OBJECT,matArray,TheMatrix,IERR) CALL MatSetFromOptions(TheMatrix,IERR) ! CALL VecCreate(PETSC_COMM_WORLD,x,IERR) CALL VecSetSizes(x,PETSC_DECIDE,2*NumberOfUnknowns,IERR) CALL VecSetFromOptions(x,IERR) CALL VecSet(x,0.+0.*PETSC_i,IERR) CALL VecSetValues(x,1,0,1.+0.*PETSC_i,INSERT_VALUES,IERR) CALL VecAssemblyBegin(x,IERR) CALL VecAssemblyEnd(x,IERR) CALL VecCreate(PETSC_COMM_WORLD,y,IERR) CALL VecSetSizes(y,PETSC_DECIDE,NumberOfUnknowns+NumberOfMeasures,IERR) CALL VecSetFromOptions(y,IERR) CALL MatMult(TheMatrix,x,y,IERR) I guess there is something wrong with the distribution of vectors between the MPI processes. With the specific sizes of my problem everything works fine with 1-2-3-6 processes, otherwise I get "Nonconforming object sizes" error, e.g. with 8 processes: [2]PETSC ERROR: - Error Message -- [2]PETSC ERROR: Nonconforming object sizes [2]PETSC ERROR: Mat mat,Vec y: local dim 383 382 [2]PETSC ERROR: Petsc Release Version 3.6.3, Dec, 03, 2015 [2]PETSC ERROR: Nonconforming object sizes ./ISS on a arch-linux2-c-debug named alessandro-HP-ProDesk-490-G2-MT by alessandro Wed Apr 6 13:31:45 2016 [2]PETSC ERROR: Configure options --with-blas-lapack-dir=/opt/intel/composer_xe_2015.2.164/mkl --with-mpi-cc=/home/alessandro/.openmpi-1.8.3/bin/mpicc --with-mpi-f90=/home/alessandro/.openmpi-1.8.3/bin/mpif90 --with-mpiexec=/home/alessandro/.openmpi-1.8.3/bin/mpiexec --with-scalar-type=complex [2]PETSC ERROR: #1 MatMult() line 2216 in /home/alessandro/CODE/PETSc/petsc-3.6.3/src/mat/interface/matrix.c ... [6]PETSC ERROR: Mat mat,Vec y: local dim 381 382 ... [7]PETSC ERROR: Mat mat,Vec y: local dim 381 382 ... [1]PETSC ERROR: Mat mat,Vec y: local dim 383 382 ... I noticed that the code works (runs and provides correct results) when the number of processes divides the length of y2 (number of rows of A21 and A22); however when I apply MatMult separately on A21 and A22 everything works fine with any number of processes Any idea where I'm messing up here? Thanks, Alessandro On Wed, 6 Apr 2016 11:16:43 +0200 Dave May wrote: On 6 April 2016 at 11:08, FRANCAVILLA MATTEO ALESSANDRO wrote: Hi, I'm trying to set a linear problem with a 2x2 block matrix (say A=[A11 A12; A21 A22] in matlab notation). Obviously I could generate a single matrix, but I would like to keep it as a 2x2 block matrix with the 4 blocks generated and used independently (e.g., because at some point I need to change some of those blocks to matrix-free). My plan was then to generate a MatShell and set a MATOP_MULT shell operation to explicitly compute y=A*x by splitting x=[x1;x2] and y=[y1;y2] and carrying out the single matrix-vector multiplications, and reassembling the final result. Take a look at MatNest - it does exactly what you want. http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatCreateNest.html http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MATNEST.html#MATNEST Thanks, Dave My problem is that my efforts in extracting x1 and x2 from the original x vector have failed. The following Fortran code is what seemed to me the most reasonable thing to do to set the matrix-vector multiplication, but it does not work !== subroutine MyMult(A,x,y,IERR) implicit none #include #include #include #include Mat A Vec x,y, x1, x2, y1, y2 IS :: ix1, ix2, iy1, iy2 PetscErrorCode ierr PetscIntM, N, II, NumberOfUnknowns, NumberOfMeasures PetscInt, ALLOCATABLE :: INDICES(:) PetscMPIInt :: RANK CALL VecGetSize(x,N,IERR) NumberOfUnknowns = N / 2 CALL VecGetSize(y,N,IERR) NumberOfMeasures = N - NumberOfUnknowns N = NumberOfUnknowns + MAX(NumberOfUnknowns,NumberOfMeasures) ALLOCATE(INDICES(N)) INDICES = (/ (ii, ii=0, N-1) /) CALL ISCreateGeneral(MPI_COMM_SELF,NumberOfUnknowns,INDICES(1:NumberOfUnknowns),PETSC_COPY_VALUES,ix1,IERR) CALL VecGetSubVector(x,ix1,x1,IERR) CALL ISCreateGeneral(MPI_COMM_SELF,NumberOfUnknowns,INDICES(NumberOfUnknowns+1:2*NumberOfUnknowns),PETSC_COPY_VALUES,ix2,IERR) CALL VecGetSubVector(x,ix2,x2,IERR) CALL ISCreateGeneral(MPI_COMM_SELF,NumberOfUnknowns,INDICES(1:NumberOfUnknowns),PETSC_COPY_VALUES,iy1,IERR) CALL VecGetSubVector(y,iy1,y1,IERR) CALL ISCreateGeneral(MPI_COMM_SELF,NumberOfUnknowns,INDICES(NumberOfUnknowns+1:NumberOfUnknowns+NumberOfMeasures),PETSC_COPY_VALUES,iy2,IERR) CALL VecGetSubVector(y,iy2,y2,IERR) CALL MatMult(L_di
[petsc-users] Correct syntax for DMDA ownership ranges in petsc4py
Hey All, Does anyone know the correct syntax for setting ownership ranges for DMDAs in petsc4py? For example, if I have a 1D DMDA with 20 vertices and 3 processes, and I want to distribute the vertices as [8,8,4], should the correct create call be: da = PETSc.DMDA().create(dim=1,dof=1,sizes=[20,],proc_sizes=[nprocs,],boundary_type=[PETSc.DM.BoundaryType.PERIODIC,], stencil_type= PETSc.DMDA.StencilType.BOX, stencil_width=1,ownership_ranges=[8,8,4]) or da = PETSc.DMDA().create(dim=1,dof=1,sizes=[20,],proc_sizes=[nprocs,],boundary_type=[PETSc.DM.BoundaryType.PERIODIC,], stencil_type= PETSc.DMDA.StencilType.BOX, stencil_width=1,ownership_ranges=[[8,8,4],]) For me, the 1st version fails with ValueError: number of dimensions 1 and number ownership ranges 3 (which I expected) but the 2nd version also fails with TypeError: Expected tuple, got list I have tried using tuples instead of lists and nothing changes. Thanks in advance for any assistance! Regards, -Chris -- Chris Eldred Postdoctoral Fellow, LAGA, University of Paris 13 PhD, Atmospheric Science, Colorado State University, 2015 DOE Computational Science Graduate Fellow (Alumni) B.S. Applied Computational Physics, Carnegie Mellon University, 2009 chris.eld...@gmail.com
Re: [petsc-users] How to extract array portions
On 6 April 2016 at 11:08, FRANCAVILLA MATTEO ALESSANDRO wrote: > Hi, > > I'm trying to set a linear problem with a 2x2 block matrix (say A=[A11 > A12; A21 A22] in matlab notation). Obviously I could generate a single > matrix, but I would like to keep it as a 2x2 block matrix with the 4 blocks > generated and used independently (e.g., because at some point I need to > change some of those blocks to matrix-free). My plan was then to generate a > MatShell and set a MATOP_MULT shell operation to explicitly compute y=A*x > by splitting x=[x1;x2] and y=[y1;y2] and carrying out the single > matrix-vector multiplications, and reassembling the final result. Take a look at MatNest - it does exactly what you want. http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatCreateNest.html http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MATNEST.html#MATNEST Thanks, Dave > My problem is that my efforts in extracting x1 and x2 from the original x > vector have failed. The following Fortran code is what seemed to me the > most reasonable thing to do to set the matrix-vector multiplication, but it > does not work > > > !== > subroutine MyMult(A,x,y,IERR) > > implicit none > > #include > #include > #include > #include > > Mat A > Vec x,y, x1, x2, y1, y2 > IS :: ix1, ix2, iy1, iy2 > PetscErrorCode ierr > PetscIntM, N, II, NumberOfUnknowns, NumberOfMeasures > PetscInt, ALLOCATABLE :: INDICES(:) > PetscMPIInt :: RANK > > CALL VecGetSize(x,N,IERR) > NumberOfUnknowns = N / 2 > CALL VecGetSize(y,N,IERR) > NumberOfMeasures = N - NumberOfUnknowns > > N = NumberOfUnknowns + MAX(NumberOfUnknowns,NumberOfMeasures) > ALLOCATE(INDICES(N)) > INDICES = (/ (ii, ii=0, N-1) /) > > CALL > ISCreateGeneral(MPI_COMM_SELF,NumberOfUnknowns,INDICES(1:NumberOfUnknowns),PETSC_COPY_VALUES,ix1,IERR) > CALL VecGetSubVector(x,ix1,x1,IERR) > CALL > ISCreateGeneral(MPI_COMM_SELF,NumberOfUnknowns,INDICES(NumberOfUnknowns+1:2*NumberOfUnknowns),PETSC_COPY_VALUES,ix2,IERR) > CALL VecGetSubVector(x,ix2,x2,IERR) > CALL > ISCreateGeneral(MPI_COMM_SELF,NumberOfUnknowns,INDICES(1:NumberOfUnknowns),PETSC_COPY_VALUES,iy1,IERR) > CALL VecGetSubVector(y,iy1,y1,IERR) > CALL > ISCreateGeneral(MPI_COMM_SELF,NumberOfUnknowns,INDICES(NumberOfUnknowns+1:NumberOfUnknowns+NumberOfMeasures),PETSC_COPY_VALUES,iy2,IERR) > CALL VecGetSubVector(y,iy2,y2,IERR) > > CALL MatMult(L_discretized,x1,y1,IERR) > CALL MatMultAdd(K_discretized,x2,y1,y1,IERR) > CALL MatMult(L_MR,x1,y2,IERR) > CALL MatMultAdd(K_MR,x2,y2,y2,IERR) > > CALL VecRestoreSubVector(y,iy1,y1,IERR) > CALL VecRestoreSubVector(y,iy2,y2,IERR) > > return > end subroutine MyMult > !== > > > Obviously the sequence of calls to ISCreateGeneral and VecGetSubVector > does not do what I expect, as the errors I'm getting are in the following > MatMult multiplications (the global sizes of x1 and x2 SHOULD BE both 771, > while y1 and y2 should be 771 and 2286): > > 1) executed with 1 MPI process: > > [0]PETSC ERROR: - Error Message > -- > [0]PETSC ERROR: Nonconforming object sizes > [0]PETSC ERROR: Mat mat,Vec y: global dim 2286 771 > [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html > for trouble shooting. > [0]PETSC ERROR: Petsc Release Version 3.6.3, Dec, 03, 2015 > [0]PETSC ERROR: ./ISS on a arch-linux2-c-debug named > alessandro-HP-ProDesk-490-G2-MT by alessandro Wed Apr 6 10:54:03 2016 > [0]PETSC ERROR: Configure options > --with-blas-lapack-dir=/opt/intel/composer_xe_2015.2.164/mkl > --with-mpi-cc=/home/alessandro/.openmpi-1.8.3/bin/mpicc > --with-mpi-f90=/home/alessandro/.openmpi-1.8.3/bin/mpif90 > --with-mpiexec=/home/alessandro/.openmpi-1.8.3/bin/mpiexec > --with-scalar-type=complex > [0]PETSC ERROR: #1 MatMult() line 2215 in > /home/alessandro/CODE/PETSc/petsc-3.6.3/src/mat/interface/matrix.c > [0]PETSC ERROR: - Error Message > -- > [0]PETSC ERROR: Nonconforming object sizes > [0]PETSC ERROR: Mat mat,Vec v3: local dim 2286 771 > [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html > for trouble shooting. > [0]PETSC ERROR: Petsc Release Version 3.6.3, Dec, 03, 2015 > [0]PETSC ERROR: ./ISS on a arch-linux2-c-debug named > alessandro-HP-ProDesk-490-G2-MT by alessandro Wed Apr 6 10:54:03 2016 > [0]PETSC ERROR: Configure options > --with-blas-lapack-dir=/opt/intel/composer_xe_2015.2.164/mkl > --with-mpi-cc=/home/alessandro/.openmpi-1.8.3/bin/mpicc > --with-mpi-f90=/home/alessandro/.openmpi-1.8.3/bin/mpif90 > --with-mpiexec=/home/alessandro/.openmpi-1.8.3/bin/mpiexec > --with-scalar-type=complex > [0]PETSC ERROR: #2 MatMultAdd() line 2396 in > /home/alessandro/CODE/PETSc/petsc-3.6.3/src/mat/interface/matrix.c > > > > 2) executed with 4 MPI process
[petsc-users] How to extract array portions
Hi, I'm trying to set a linear problem with a 2x2 block matrix (say A=[A11 A12; A21 A22] in matlab notation). Obviously I could generate a single matrix, but I would like to keep it as a 2x2 block matrix with the 4 blocks generated and used independently (e.g., because at some point I need to change some of those blocks to matrix-free). My plan was then to generate a MatShell and set a MATOP_MULT shell operation to explicitly compute y=A*x by splitting x=[x1;x2] and y=[y1;y2] and carrying out the single matrix-vector multiplications, and reassembling the final result. My problem is that my efforts in extracting x1 and x2 from the original x vector have failed. The following Fortran code is what seemed to me the most reasonable thing to do to set the matrix-vector multiplication, but it does not work !== subroutine MyMult(A,x,y,IERR) implicit none #include #include #include #include Mat A Vec x,y, x1, x2, y1, y2 IS :: ix1, ix2, iy1, iy2 PetscErrorCode ierr PetscIntM, N, II, NumberOfUnknowns, NumberOfMeasures PetscInt, ALLOCATABLE :: INDICES(:) PetscMPIInt :: RANK CALL VecGetSize(x,N,IERR) NumberOfUnknowns = N / 2 CALL VecGetSize(y,N,IERR) NumberOfMeasures = N - NumberOfUnknowns N = NumberOfUnknowns + MAX(NumberOfUnknowns,NumberOfMeasures) ALLOCATE(INDICES(N)) INDICES = (/ (ii, ii=0, N-1) /) CALL ISCreateGeneral(MPI_COMM_SELF,NumberOfUnknowns,INDICES(1:NumberOfUnknowns),PETSC_COPY_VALUES,ix1,IERR) CALL VecGetSubVector(x,ix1,x1,IERR) CALL ISCreateGeneral(MPI_COMM_SELF,NumberOfUnknowns,INDICES(NumberOfUnknowns+1:2*NumberOfUnknowns),PETSC_COPY_VALUES,ix2,IERR) CALL VecGetSubVector(x,ix2,x2,IERR) CALL ISCreateGeneral(MPI_COMM_SELF,NumberOfUnknowns,INDICES(1:NumberOfUnknowns),PETSC_COPY_VALUES,iy1,IERR) CALL VecGetSubVector(y,iy1,y1,IERR) CALL ISCreateGeneral(MPI_COMM_SELF,NumberOfUnknowns,INDICES(NumberOfUnknowns+1:NumberOfUnknowns+NumberOfMeasures),PETSC_COPY_VALUES,iy2,IERR) CALL VecGetSubVector(y,iy2,y2,IERR) CALL MatMult(L_discretized,x1,y1,IERR) CALL MatMultAdd(K_discretized,x2,y1,y1,IERR) CALL MatMult(L_MR,x1,y2,IERR) CALL MatMultAdd(K_MR,x2,y2,y2,IERR) CALL VecRestoreSubVector(y,iy1,y1,IERR) CALL VecRestoreSubVector(y,iy2,y2,IERR) return end subroutine MyMult !== Obviously the sequence of calls to ISCreateGeneral and VecGetSubVector does not do what I expect, as the errors I'm getting are in the following MatMult multiplications (the global sizes of x1 and x2 SHOULD BE both 771, while y1 and y2 should be 771 and 2286): 1) executed with 1 MPI process: [0]PETSC ERROR: - Error Message -- [0]PETSC ERROR: Nonconforming object sizes [0]PETSC ERROR: Mat mat,Vec y: global dim 2286 771 [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting. [0]PETSC ERROR: Petsc Release Version 3.6.3, Dec, 03, 2015 [0]PETSC ERROR: ./ISS on a arch-linux2-c-debug named alessandro-HP-ProDesk-490-G2-MT by alessandro Wed Apr 6 10:54:03 2016 [0]PETSC ERROR: Configure options --with-blas-lapack-dir=/opt/intel/composer_xe_2015.2.164/mkl --with-mpi-cc=/home/alessandro/.openmpi-1.8.3/bin/mpicc --with-mpi-f90=/home/alessandro/.openmpi-1.8.3/bin/mpif90 --with-mpiexec=/home/alessandro/.openmpi-1.8.3/bin/mpiexec --with-scalar-type=complex [0]PETSC ERROR: #1 MatMult() line 2215 in /home/alessandro/CODE/PETSc/petsc-3.6.3/src/mat/interface/matrix.c [0]PETSC ERROR: - Error Message -- [0]PETSC ERROR: Nonconforming object sizes [0]PETSC ERROR: Mat mat,Vec v3: local dim 2286 771 [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting. [0]PETSC ERROR: Petsc Release Version 3.6.3, Dec, 03, 2015 [0]PETSC ERROR: ./ISS on a arch-linux2-c-debug named alessandro-HP-ProDesk-490-G2-MT by alessandro Wed Apr 6 10:54:03 2016 [0]PETSC ERROR: Configure options --with-blas-lapack-dir=/opt/intel/composer_xe_2015.2.164/mkl --with-mpi-cc=/home/alessandro/.openmpi-1.8.3/bin/mpicc --with-mpi-f90=/home/alessandro/.openmpi-1.8.3/bin/mpif90 --with-mpiexec=/home/alessandro/.openmpi-1.8.3/bin/mpiexec --with-scalar-type=complex [0]PETSC ERROR: #2 MatMultAdd() line 2396 in /home/alessandro/CODE/PETSc/petsc-3.6.3/src/mat/interface/matrix.c 2) executed with 4 MPI processes (I'm copying the output of process 0 only, they're all the same): [0]PETSC ERROR: - Error Message -- [0]PETSC ERROR: Nonconforming object sizes [0]PETSC ERROR: Mat mat,Vec x: global dim 771 3084 [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting. [0]PETSC ERROR: Petsc Release Version 3.6.3, Dec, 03, 2015 [0]PETSC ERROR: ./ISS on a arc
