using MatCreateSeqBDiag

[Alejandro Garzon]

Hi, I'm trying to use MatCreateSeqBDiag. I wrote a small test program
to see how this function works. The program creates a matrix A,
prints it in matlab ascii format and in binary format, then creates a
vector x, calculates A*x and puts the result in vector b.

I have found the following inconsistencies:
1) While I intended the matrix to have real elements, it is printed in
matlab ascii format as if all the elements were imaginary. The
elements are correct except for an extra "i" to the right of all of
them.

2) When printed in binary format, the matrix appears with all its
elements set to zero.

3) The internal representation of the matrix seems to be Ok as a
matrix-vector multiplication A*x gives the right result (without the
i's).

The source code and the output of the program are below. Can you take
a look at it and tell me if the I haven't built the matrix in the
right way or if it could be a bug in the MatView function when
handling this matrix format ? Thanks.

// ************** Source code program *********************
// This code is for checking the usage of MatCreateSeqBDiag

#include<stdio.h>
#include<stdlib.h>
#include"petscksp.h"

int main(int argc,char **args)
{
  Mat A;
  Vec x,b;
  PetscInt diag[9],ixm[10];
  PetscScalar **diagv,*diagv1,vval[10];
  int i;
  PetscViewer v;
  PetscErrorCode ierr;

  PetscInitialize(&argc,&args,(char *)0,(char *)0);

  // A 10x10 matrix with 2x2 blocks will be build.
  // Only blocks in the main diagonal will be filled.
  diag[0]=0;
  diagv=malloc(1*sizeof(PetscScalar*));
  diagv1=malloc(20*sizeof(PetscScalar));

  // The main block diagonal will contain the numbers from
  // 1 to 20
  for(i=0;i<20;i++){
    diagv1[i]=(i+1.0);
  }

  diagv[0]=diagv1;

  ierr=MatCreateSeqBDiag(PETSC_COMM_WORLD,10,10,1,2,diag,diagv,&A);
  CHKERRQ(ierr);
  ierr=MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr=MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

  // Printing matrix A to the stardard output in matlab ascii format.
  ierr=PetscViewerSetFormat(PETSC_VIEWER_STDOUT_WORLD,
			    PETSC_VIEWER_ASCII_MATLAB);CHKERRQ(ierr);
  ierr=MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

  // Printing matrix A to binary file "binary_file"
  ierr=PetscViewerBinaryOpen(PETSC_COMM_WORLD,"binary_file",
			     FILE_MODE_WRITE,&v);CHKERRQ(ierr);
  ierr=MatView(A,v);CHKERRQ(ierr);

  // Building vector with all elements equal to 1
  ierr=VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);
  ierr=VecSetSizes(x,PETSC_DECIDE,10);CHKERRQ(ierr);
  ierr=VecSetFromOptions(x);CHKERRQ(ierr);

  for(i=0;i<10;i++){
    ixm[i]=i;
    vval[i]=1;
  }

  ierr=VecSetValues(x,10,ixm,vval,INSERT_VALUES);CHKERRQ(ierr);
  ierr=VecAssemblyBegin(x);CHKERRQ(ierr);
  ierr=VecAssemblyEnd(x);CHKERRQ(ierr);

  // Printing vector x to standard output in matlab ascii format
  ierr=VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

  ierr=VecDuplicate(x,&b);CHKERRQ(ierr);

  // b = A*x
  ierr=MatMult(A,x,b);CHKERRQ(ierr);

  // Printing vector b to standard output in matlab ascii format
  ierr=VecView(b,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

  PetscFinalize();

  return 0;
}

// ****************** Program output *********************
% Size = 10 10
% Nonzeros = 20
zzz = zeros(20,3);
zzz = [
1 1  1.0000000000000000e+00i
1 2  3.0000000000000000e+00i
2 1  2.0000000000000000e+00i
2 2  4.0000000000000000e+00i
3 3  5.0000000000000000e+00i
3 4  7.0000000000000000e+00i
4 3  6.0000000000000000e+00i
4 4  8.0000000000000000e+00i
5 5  9.0000000000000000e+00i
5 6  1.1000000000000000e+01i
6 5  1.0000000000000000e+01i
6 6  1.2000000000000000e+01i
7 7  1.3000000000000000e+01i
7 8  1.5000000000000000e+01i
8 7  1.4000000000000000e+01i
8 8  1.6000000000000000e+01i
9 9  1.7000000000000000e+01i
9 10  1.9000000000000000e+01i
10 9  1.8000000000000000e+01i
10 10  2.0000000000000000e+01i
];
 Mat_0 = spconvert(zzz);
Vec_1 = [
1.0000000000000000e+00
1.0000000000000000e+00
1.0000000000000000e+00
1.0000000000000000e+00
1.0000000000000000e+00
1.0000000000000000e+00
1.0000000000000000e+00
1.0000000000000000e+00
1.0000000000000000e+00
1.0000000000000000e+00
];
Vec_2 = [
4.0000000000000000e+00
6.0000000000000000e+00
1.2000000000000000e+01
1.4000000000000000e+01
2.0000000000000000e+01
2.2000000000000000e+01
2.8000000000000000e+01
3.0000000000000000e+01
3.6000000000000000e+01
3.8000000000000000e+01
];

--
Alejandro