[petsc-users] MATELEMENTAL strange behaviour

[Bernard Knaepen]

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/7f3c2d89/attachment.html>