[petsc-users] MATELEMENTAL strange behaviour

[Hong]

Bernard:
With your code, I reproduced error and found that the default
MAT_ROW_ORIENTED for MatSetValues() is changed to MAT_COLUMN_ORIENTED. This
is a bug in our library. I'll fix it.

You can add
MatSetOption(C,MAT_ROW_ORIENTED,PETSC_TRUE);
before 2nd call of matrixC(C,t).

Thanks for reporting the bug.

Hong

Hello,
>
> Here is a trimmed down piece of code that illustrates a strange behaviour
> that I am desperately trying to understand.
>
> The problem is that I would expect successive calls to MatView(C,…) to
> display the same numbers. However, the numbers displayed are different,
> even though the content of the matrix C should be the same.
>
> Any help to resolve this would be appreciated.
>
> Cheers,
> Bernard.
>
> Code:
> ******
>
> static char help[] = "Elemental test\n\n";
>
> /*T
>
> T*/
>
> #include <petscmat.h>
>
> PetscScalar pi=3.141592653589793;
>
> /* Prototypes */
> PetscErrorCode matrixC(Mat C, PetscScalar t);
>
> #undef __FUNCT__
> #define __FUNCT__ "main"
> int main(int argc, char **args){
>
>    Mat C;
>
>    PetscInt N = 8;
>    PetscScalar t=0.1, dt=0.0001;
>    PetscErrorCode ierr;
>
>    /* parameters */
>    ierr = PetscInitialize(&argc, &args, (char*) 0, help) ;CHKERRQ(ierr);
>
>    ierr = MatCreate(PETSC_COMM_WORLD, &C); CHKERRQ(ierr);
>    ierr = MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, N, N); CHKERRQ(ierr);
>    ierr = MatSetType(C, MATELEMENTAL); CHKERRQ(ierr);
>    ierr = MatSetFromOptions(C); CHKERRQ(ierr);
>    ierr = MatSetUp(C); CHKERRQ(ierr);
>
>    ierr=matrixC(C,t);
>    ierr = MatView(C,PETSC_VIEWER_STDOUT_WORLD);
>
>    ierr=matrixC(C,t);
>    ierr = MatView(C,PETSC_VIEWER_STDOUT_WORLD);
>
>    ierr = MatDestroy(&C);CHKERRQ(ierr);
>
>    ierr = PetscFinalize();
>
>    return ierr;
> }
>
> /* Matrix C*/
> PetscErrorCode matrixC(Mat C, PetscScalar t){
>
>    PetscErrorCode ierr;
>    IS isrows,iscols;
>    const PetscInt *rows,*cols;
>    PetscScalar    *v;
>    PetscInt i,j,nrows,ncols;
>    PetscInt n,m;
>
>    /* Set local matrix entries */
>    ierr = MatGetOwnershipIS(C,&isrows,&iscols);CHKERRQ(ierr);
>    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++) {
>       n=rows[i];
>       for (j=0; j<ncols; j++) {
>          m=cols[j];
>          v[i*ncols+j] = -0.5*(exp(-(m+n+1.5)*(m+n+1.5)
> *pi*pi*t)-exp(-(m-n+0.5)*(m-n+0.5)*pi*pi*t));
>       }
>    }
>
>    ierr = MatSetValues(C,nrows,rows,ncols,cols,v,INSERT_VALUES);
> CHKERRQ(ierr);
>    ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
>    ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
>    ierr = ISRestoreIndices(isrows,&rows);CHKERRQ(ierr);
>    ierr = ISRestoreIndices(iscols,&cols);CHKERRQ(ierr);
>    ierr = ISDestroy(&isrows);CHKERRQ(ierr);
>    ierr = ISDestroy(&iscols);CHKERRQ(ierr);
>    ierr = PetscFree(v);CHKERRQ(ierr);
>
>     PetscFunctionReturn(0);
> }
>
> Output:
> ********
>
> mpirun -n 2 ./elemental
>
> Mat Object: 2 MPI processes
>   type: elemental
> Elemental matrix (cyclic ordering)
> 0.336403 2.8059e-06 0.0532215 1.04511e-09 0.00104438
> 0.00104718 0.390672 0.0542687 0.0542687 0.390672
> 0.389625 0.00104718 0.390669 2.80695e-06 0.0542687
> 2.80695e-06 0.390672 0.00104718 0.390672 0.0542687
> 0.0542659 0.0542687 0.390672 0.00104718 0.390672
>
> Elemental matrix (explicit ordering)
>   Mat Object: 2 MPI processes
>     type: mpidense
> 3.3640319378153799e-01 5.3221486768492490e-02 1.0443777750910219e-03
> 2.8059034408741489e-06 1.0451056906250189e-09
> 3.8962468055003047e-01 3.9066905832512150e-01 5.4268670447024388e-02
> 1.0471847236375868e-03 2.8069486006233780e-06
> 5.4265864543583515e-02 3.9067186422856237e-01 3.9067186527366804e-01
> 5.4268671492184138e-02 1.0471847236916457e-03
> 1.0471836785318962e-03 5.4268671492130077e-02 3.9067186527372211e-01
> 3.9067186527372211e-01 5.4268671492184138e-02
> 2.8069485465647739e-06 1.0471847236916453e-03 5.4268671492184138e-02
> 3.9067186527372211e-01 3.9067186527372211e-01
> Mat Object: 2 MPI processes
>   type: elemental
> Elemental matrix (cyclic ordering)
> 0.336403 2.8059e-06 0.390672 2.80695e-06 0.00104718
> 0.0532215 1.04511e-09 0.389625 0.390672 0.0542687
> 0.00104438 0.390672 0.390669 0.390672 0.0542659
> 0.00104718 0.0542687 0.0542687 0.0542687 0.390672
> 0.0542687 0.00104718 2.80695e-06 0.00104718 0.390672
>
> Elemental matrix (explicit ordering)
>   Mat Object: 2 MPI processes
>     type: mpidense
> 3.3640319378153799e-01 3.9067186527372211e-01 1.0471847236916453e-03
> 2.8059034408741489e-06 2.8069486006233780e-06
> 1.0443777750910219e-03 3.9066905832512150e-01 5.4265864543583515e-02
> 3.9067186527372211e-01 3.9067186527372211e-01
> 5.4268671492130077e-02 2.8069485465647739e-06 3.9067186527366804e-01
> 1.0471847236375868e-03 1.0471847236916457e-03
> 5.3221486768492490e-02 3.8962468055003047e-01 5.4268671492184138e-02
> 1.0451056906250189e-09 3.9067186527372211e-01
> 1.0471836785318962e-03 5.4268670447024388e-02 3.9067186422856237e-01
> 5.4268671492184138e-02 5.4268671492184138e-02
>
>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20160928/02d35a83/attachment-0001.html>