using MatCreateSeqBDiag

[Matthew Knepley]

On 8/28/07, Alejandro Garzon <gtg100n at mail.gatech.edu> wrote:
> 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.

This was a cut&paste error. Fixed.

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

Fixed this now too.

These fixes will appear in the next patch release. Hopefully before the end of
the week.

  Thanks for reporting these,

    Matt

> 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
>
>
>
>
>


-- 
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