[petsc-users] MATMPIBAIJ matrix allocation

[Jose E. Roman]

El 29/07/2013, a las 14:54, Heikki Virtanen escribió:

> Hi, I try to solve an eigenvalue problem (Ax = lambda B x) where A 
> and B are complex matrices. Unfortunately, my data structures 
> are not capable of handling directly complex numbers, yet. So,
> I have to use 
> 
>                 [re(A) -im(A)]
> RealA  = [                    ]
>                 [imA)   re(A)] 
> 
> matrices instead. How should I allocate matrices if I know
> compressed sparse row matrix formats of A and B matrices. I have used
> something like this.
> 
> ierr = MatCreate(PETSC_COMM_WORLD,&RealA); CHKERRQ(ierr);    
> ierr = MatSetSizes(RealA,PETSC_DECIDE,PETSC_DECIDE,2*n,2*n); CHKERRQ(ierr);
> ierr = MatSetType(RealA,MATMPIBAIJ); CHKERRQ(ierr);
> ierr = MatSetBlockSize(RealA,2);CHKERRQ(ierr);
> ierr = MatSetFromOptions(RealA);CHKERRQ(ierr);
>                 
> ierr = MatMPIBAIJSetPreallocationCSR (RealA,2,rows,cols,0); CHKERRQ(ierr);
> ierr = MatAssemblyBegin (RealA,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
> ierr = MatAssemblyEnd (RealA,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
> 
> where n is the global size of matrix A. rows and cols are csr arrays of A.
> Each submatrix of RealA matrix have the same nonzero pattern than A.
> But I am not sure if this is correct, because when I print my matrices out
> they look more like band matrices than block matrices.
> 
> -Heikki

You are creating a BAIJ matrix with the nonzero pattern defined by (rows,cols), where each entry of the matrix is a 2x2 block, in your case
           [ re(a_ij) -im(a_ij) ]
           [ im(a_ij)  re(a_ij) ]

So the matrix you are creating is not this one
                [re(A) -im(A)]
                [imA)   re(A)]
but the one resulting from a perfect shuffle permutation. That's why you are not seeing a 2x2 block structure for the whole matrix.

Jose