[petsc-users] MATELEMENTAL MatSetValue( )
Jed Brown
jed at jedbrown.org
Thu May 26 22:54:33 CDT 2016
Please always use "reply-all" so that your messages go to the list.
This is standard mailing list etiquette. It is important to preserve
threading for people who find this discussion later and so that we do
not waste our time re-answering the same questions that have already
been answered in private side-conversations. You'll likely get an
answer faster that way too.
"Eric D. Robertson" <eric.robertson at cfdrc.com> writes:
> Hi Jed,
>
> I actually tried your suggestion (which is in-line with ex38.c). Running that example actually seems to produce similar behavior. The 'C' matrix actually seems to grow in size with the number of processors. Is there something wrong here, or am I misinterpreting the purpose of the example?
That test (NOT tutorial) uses
ierr = MatSetSizes(C,m,n,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
which keeps the local sizes constant. You should be creating the matrix
size that is semantically meaningful to you and setting the owned rows
as identified by MatGetOwnershipIS().
> Eric Robertson
> Research Engineer, BET
> CFD Research Corporation
> 701 McMillian Way NW, Suite D
> Huntsville, AL 35806
> Ph: 256-726-4912
> Email: eric.robertson at cfdrc.com
>
>
> ________________________________________
> From: Jed Brown [jed at jedbrown.org]
> Sent: Wednesday, May 25, 2016 10:32 PM
> To: Eric D. Robertson; petsc-users at mcs.anl.gov
> Subject: Re: [petsc-users] MATELEMENTAL MatSetValue( )
>
> Use MatGetOwnershipIS() and set only the entries in your owned rows and
> columns. What you're doing is nonscalable anyway because every process
> is filling the global matrix. If you want a cheap hack for small
> problems on small numbers of cores, you can run your insertion loop on
> only rank 0.
>
> "Eric D. Robertson" <eric.robertson at cfdrc.com> writes:
>
>> I am trying to multiply two dense matrices using the Elemental interface. I fill the matrix using MatSetValue( ) like below:
>>
>> for ( i = 0; i < Matrix.M; i++){
>> for ( j = 0; j < Matrix.N; j++) {
>> PetscScalar temp = i + one + (j*three);
>> MatSetValue(Matrix.A, i, j, temp, INSERT_VALUES);
>>
>> }
>> }
>>
>>
>> However, I seem to get the following error:
>>
>> [0]PETSC ERROR: No support for this operation for this object type
>> [0]PETSC ERROR: Only ADD_VALUES to off-processor entry is supported
>>
>> But if I use ADD_VALUES, I get a different matrix depending on the number of processors used. The entries become multiplied by the number of processors. How do I reconcile this?
>>
>>
>> <redir.aspx?REF=hVEx-ekIltbV6NDJuYjzplTZG3ymEirNPcO9eR2dIGrvJzCyE4XTCAFtYWlsdG86ZXJpYy5yb2JlcnRzb25AY2ZkcmMuY29t>
>
> static char help[] = "Test interface of Elemental. \n\n";
>
> #include <petscmat.h>
>
> #undef __FUNCT__
> #define __FUNCT__ "main"
> int main(int argc,char **args)
> {
> Mat C,Caij;
> PetscInt i,j,m = 5,n,nrows,ncols;
> const PetscInt *rows,*cols;
> IS isrows,iscols;
> PetscErrorCode ierr;
> PetscBool flg,Test_MatMatMult=PETSC_FALSE,mats_view=PETSC_FALSE;
> PetscScalar *v;
> PetscMPIInt rank,size;
>
> PetscInitialize(&argc,&args,(char*)0,help);
> ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
> ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
> ierr = PetscOptionsHasName(NULL,"-mats_view",&mats_view);CHKERRQ(ierr);
>
> /* Get local block or element size*/
> ierr = PetscOptionsGetInt(NULL,"-m",&m,NULL);CHKERRQ(ierr);
> n = m;
> ierr = PetscOptionsGetInt(NULL,"-n",&n,NULL);CHKERRQ(ierr);
>
> ierr = MatCreate(PETSC_COMM_WORLD,&C);CHKERRQ(ierr);
> ierr = MatSetSizes(C,m,n,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
> ierr = MatSetType(C,MATELEMENTAL);CHKERRQ(ierr);
> ierr = MatSetFromOptions(C);CHKERRQ(ierr);
> ierr = MatSetUp(C);CHKERRQ(ierr);
>
> ierr = PetscOptionsHasName(NULL,"-row_oriented",&flg);CHKERRQ(ierr);
> if (flg) {ierr = MatSetOption(C,MAT_ROW_ORIENTED,PETSC_TRUE);CHKERRQ(ierr);}
> ierr = MatGetOwnershipIS(C,&isrows,&iscols);CHKERRQ(ierr);
> ierr = PetscOptionsHasName(NULL,"-Cexp_view_ownership",&flg);CHKERRQ(ierr);
> if (flg) { /* View ownership of explicit C */
> IS tmp;
> ierr = PetscPrintf(PETSC_COMM_WORLD,"Ownership of explicit C:\n");CHKERRQ(ierr);
> ierr = PetscPrintf(PETSC_COMM_WORLD,"Row index set:\n");CHKERRQ(ierr);
> ierr = ISOnComm(isrows,PETSC_COMM_WORLD,PETSC_USE_POINTER,&tmp);CHKERRQ(ierr);
> ierr = ISView(tmp,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
> ierr = ISDestroy(&tmp);CHKERRQ(ierr);
> ierr = PetscPrintf(PETSC_COMM_WORLD,"Column index set:\n");CHKERRQ(ierr);
> ierr = ISOnComm(iscols,PETSC_COMM_WORLD,PETSC_USE_POINTER,&tmp);CHKERRQ(ierr);
> ierr = ISView(tmp,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
> ierr = ISDestroy(&tmp);CHKERRQ(ierr);
> }
>
> /* Set local matrix entries */
> ierr = ISGetLocalSize(isrows,&nrows);CHKERRQ(ierr);
> ierr = ISGetIndices(isrows,&rows);CHKERRQ(ierr);
> ierr = ISGetLocalSize(iscols,&ncols);CHKERRQ(ierr);
> ierr = ISGetIndices(iscols,&cols);CHKERRQ(ierr);
> ierr = PetscMalloc1(nrows*ncols,&v);CHKERRQ(ierr);
> for (i=0; i<nrows; i++) {
> for (j=0; j<ncols; j++) {
> /*v[i*ncols+j] = (PetscReal)(rank);*/
> v[i*ncols+j] = (PetscReal)(rank*10000+100*rows[i]+cols[j]);
> }
> }
> ierr = MatSetValues(C,nrows,rows,ncols,cols,v,INSERT_VALUES);CHKERRQ(ierr);
> ierr = ISRestoreIndices(isrows,&rows);CHKERRQ(ierr);
> ierr = ISRestoreIndices(iscols,&cols);CHKERRQ(ierr);
> ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
> ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
>
> /* Test MatView() */
> if (mats_view) {
> ierr = MatView(C,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
> }
>
> /* Set unowned matrix entries - add subdiagonals and diagonals from proc[0] */
> if (!rank) {
> PetscInt M,N,cols[2];
> ierr = MatGetSize(C,&M,&N);CHKERRQ(ierr);
> for (i=0; i<M; i++) {
> cols[0] = i; v[0] = i + 0.5;
> cols[1] = i-1; v[1] = 0.5;
> if (i) {
> ierr = MatSetValues(C,1,&i,2,cols,v,ADD_VALUES);CHKERRQ(ierr);
> } else {
> ierr = MatSetValues(C,1,&i,1,&i,v,ADD_VALUES);CHKERRQ(ierr);
> }
> }
> }
> ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
> ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
>
> /* Test MatMult() */
> ierr = MatComputeExplicitOperator(C,&Caij);CHKERRQ(ierr);
> ierr = MatMultEqual(C,Caij,5,&flg);CHKERRQ(ierr);
> if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"C != Caij. MatMultEqual() fails");
> ierr = MatMultTransposeEqual(C,Caij,5,&flg);CHKERRQ(ierr);
> if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"C != Caij. MatMultTransposeEqual() fails");
>
> /* Test MatMultAdd() and MatMultTransposeAddEqual() */
> ierr = MatMultAddEqual(C,Caij,5,&flg);CHKERRQ(ierr);
> if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"C != Caij. MatMultAddEqual() fails");
> ierr = MatMultTransposeAddEqual(C,Caij,5,&flg);CHKERRQ(ierr);
> if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"C != Caij. MatMultTransposeAddEqual() fails");
>
> /* Test MatMatMult() */
> ierr = PetscOptionsHasName(NULL,"-test_matmatmult",&Test_MatMatMult);CHKERRQ(ierr);
> if (Test_MatMatMult) {
> Mat CCelem,CCaij;
> ierr = MatMatMult(C,C,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&CCelem);CHKERRQ(ierr);
> ierr = MatMatMult(Caij,Caij,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&CCaij);CHKERRQ(ierr);
> ierr = MatMultEqual(CCelem,CCaij,5,&flg);CHKERRQ(ierr);
> if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"CCelem != CCaij. MatMatMult() fails");
> ierr = MatDestroy(&CCaij);CHKERRQ(ierr);
> ierr = MatDestroy(&CCelem);CHKERRQ(ierr);
> }
>
> ierr = PetscFree(v);CHKERRQ(ierr);
> ierr = ISDestroy(&isrows);CHKERRQ(ierr);
> ierr = ISDestroy(&iscols);CHKERRQ(ierr);
> ierr = MatDestroy(&C);CHKERRQ(ierr);
> ierr = MatDestroy(&Caij);CHKERRQ(ierr);
> ierr = PetscFinalize();
> return 0;
> }
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 818 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20160526/a0e911d0/attachment.pgp>
More information about the petsc-users
mailing list