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
